stat.AP 2015-03-14 0

A bivariate space-time downscaler under space and time misalignment

Ozone and particulate matter PM2.5 are co-pollutants that have long been associated with increased public health risks. Information on concentration levels for both pollutants come from two sources: monitoring sites and output from complex numerical …

Authors: Veronica J. Berrocal, Alan E. Gelf, David M. Holl

Submitte d to the Annals of Applie d Statistics A BIV ARIA TE SP A CE-TIME DO WNSCALER UNDER SP A CE AND TIME MISALIGNMENT By Vero nica J. Berrocal ∗ , Alan E. Gelf and † and D a vid M. Holland ‡ SAMSI ∗ , Duke University † and U.S. Envir onmental Pr ote ction A gency ‡ Ozone and particulate matter PM 2 . 5 are co-p ollutan ts that hav e long been asso ciated with increased public health risks. Information on concentration lev els for b oth p ollutan ts come from t w o sources: monitoring sites and output from complex numerical models that pro- duce concentration surfaces o v er large spatial regions. In this pap er, w e offer a fully-model based approach for fusing these t wo sources of information for the pair of co-p ollutan ts whic h is computationally feasible ov er large spatial regions and long p eriods of time. Due to the asso ciation betw een concentration levels of the t wo environmen- tal contaminan ts, it is exp ected that information regarding one will help to improv e prediction of the other. Misalignment is an ob vious issue since the monitoring net w orks for the tw o contaminan ts only partly intersect and because the collection rate for PM 2 . 5 is typically less frequent than that for ozone. Extending previous work in Berro cal et al. (2010), we introduce a biv ariate do wnscaler that provides a flexible class of biv ariate space- time assimilation mo dels. W e discuss computational issues for mo del fitting and analyze a dataset for ozone and PM 2 . 5 for the ozone sea- son during year 2002. W e sho w a mo dest impro vemen t in predictive p erformance, not surprising in a setting where w e can an ticipate only a small gain. 1. In tro duction. Ozone and particulate matter, PM 2 . 5 , hav e long b een asso ciated with increased public health risks, e.g., of respiratory diseases [ 4 , 12 , 45 ], cardio v ascular diseases [ 4 , 12 ], and mortality and morbidity in general [ 13 , 47 ]. T o set air quality standards the EP A utilizes information from monitoring net w orks and from air quality computer mo dels. The tw o sources of information are v aluable in differen t wa ys. The ob- serv ations rep orted b y monitoring netw orks, though sparsely collected and, sometimes, affected b y missingness, provide direct measuremen ts of the true v alue up to measurement error. The output from air quality mo dels esti- mates p ollutan t concentrations ov er large spatial domains at the grid cell AMS 2000 subje ct classifications: Primary 60K35, 60K35; secondary 60K35 Keywor ds and phrases : Co-kriging, coregionalization, dynamic mo del, kriging, m ulti- v ariate spatial pro cess, spatially v arying co efficien ts 1 imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 2 V. J. BERR OCAL ET AL. lev el. The estimates are view ed as a verages o ver these cells and do not con- tain any missing data, but are uncalibrated. It is desirable to com bine both sources of information but to do so, it is necessary to solv e the inherent “c hange of supp ort” problem of the data [ 2 , 9 , 22 ], that is, the spatial mis- alignmen t b et ween the observ ational data and the n umerical mo del output. Addressing the difference in spatial resolution b et w een the tw o sources of data will also allow to ev aluate and calibrate the numerical mo del output. W e are not in terested in calibration of the n umerical mo del in the sense of [ 29 ]. Rather our fo cus is on resolving the difference in spatial scale be- t w een the numerical mo del output and the observ ational data, while bias- correcting the predictions generated by a n umerical mo del. In this regard, sev eral metho ds ha v e been prop osed to assess n umerical mo dels and address the issue of “incommensurabilit y” b et w een grid mo del av erages and p oin t measuremen ts [ 50 ]. Such effort is also common in the con text of data as- similation within the atmospheric science literature, where the goal is to com bine observ ational data on the current state of the atmosphere with a short-range forecast produced by a n umerical w eather prediction model to yield initial conditions for a n umerical atmospheric mo del. The approaches tend to b e algorithmic and ad ho c, only o ccasionally model-based, and do not necessarily address do wnscaling. F uller discussion can b e found in [ 10 ] and [ 28 ]. In the context of air quality , [ 36 ] prop ose to ev aluate the hourly ozone predictions generated b y a geophysical mo del on the grid-cell scale b y using the observ ations tak en at monitoring sites, and predicting hourly lev el of ozone at the mo del grid cells by a v eraging the predictions at M regularly-spaced sites taken within eac h grid cell (see, as well, [ 15 ]). A dif- feren t strategy has been prop osed by [ 27 ] who ignore the difference in spatial resolution b et w een mo del output and observ ations, rather suggesting to ev al- uate a numerical mo del b y lo oking for differences b et w een the mo del output and the observ ations in terms of v ariograms and correlograms. [ 16 ] dev elop a “fusion” model, expressing the n umerical mo del output as the integral ov er a grid cell (scaled by the area of the cell) of a latent p oin t- lev el pro cess. Their main goal is to recov er the true unobserved pro cess, but, as a by-product, they obtain estimates of the bias of the numerical mo del output, thus allowing ev aluation and calibration of the numerical mo del. The work is an instance of Ba y esian melding [ 39 ], and has gained considerable p opularit y [ 14 , 46 ]. See also recent work in this spirit applied to sulfate aerosol [ 49 ] and ammonium [ 11 ]. The Bay esian melding approac h of [ 16 ], though p opular, suffers from sev eral limitations, as p oin ted out b y [ 33 ] and [ 3 ]. Additionally , it is only computationally feasible as a spatial mo del. A spatio-temporal extension, prop osed by [ 35 ], sp ecifies the latent imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 3 pro cess at the grid cell lev el, leading to comparison b et w een the n umerical mo del output and the observ ational data at the grid cell scale. Op erating at the p oin t lev el, [ 23 ] propose to correct the mo del output by do wnscaling the predictions from grid cells to p oin ts and comparing them with the observ ations. F or each site, a tw o-step linear regression mo del that relates the observ ed hourly ozone level at a site with the n umerical mo del output at the grid cell in which the site lies is prop osed. Building up on the work of [ 23 ], [ 32 ] suggest interpolating the in tercept and the co efficien t of the numerical mo del output, estimated at each site, via kriging, thus dev eloping a “spatio-temp oral” mo del that will allo w c orr e cte d prediction of ozone lev el also at unmonitored sites. Ho w ever, this ad-ho c approac h does not provide estimates of the uncertainties asso ciated with the predictions. A formal spatio-temp oral model, that builds up on the do wnscaling idea of [ 23 ], has b een prop osed as a univ ariate do wnscaler b y [ 3 ] with a version implemen ted and describ ed in [ 42 ] and in [ 18 ]. The con tribution of this paper is to extend the spatio-temp oral downscaler from the univ ariate setting to a biv ariate setting in order to exploit the cor- relation b oth in the observed concen tration lev els of ozone and PM 2 . 5 and in the concen tration lev els of ozone and PM 2 . 5 , pro vided b y the n umerical mo del. W orking in a biv ariate setting will prov e particularly adv antageous for predicting particulate matter. The sampling frequency of the PM 2 . 5 mon- itors - most of the monitors measure PM 2 . 5 concen tration ev ery 3 da ys - yields c hallenging in terp olation of the monitoring data to the entire spatial domain as w ell as difficult ev aluation of the predictions generated b y the n u- merical mo del. Rather, in en vironmental health studies, researc hers alw ays use monitoring data to characterize exp osure to fine particulate matter, dealing with the PM 2 . 5 sampling frequency issue by aggregating the mon- itoring data ov er time and/or ov er space. Similarly , most of the research effort on PM 2 . 5 has b een devoted to dev eloping spatial and spatio-temp oral mo dels to predict PM 2 . 5 based on monitoring data, meteorological co v ari- ates or observed concentrations for another p ollutan t [ 5 , 30 , 31 , 40 , 43 , 48 ]. A fusion mo del to combine monitoring data for PM 2 . 5 with satellite aerosol optical depth (AOD) data has b een prop osed by [ 37 ]: ho w ever, the analysis w as conducted on monthly av erage concentration rather than daily and it rev ealed that A OD pro vides no impro v ement in predicting fine particulate matter. In this paper, w e presen t a general biv ariate spatio-temporal do wnscaler dev elop ed for the numerical mo del outputs of an air quality mo del; we illus- trate the metho dology with regard to prediction of concen tration of ozone and fine particulate matter. The mo del, which to our knowledge, is the imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 4 V. J. BERR OCAL ET AL. first biv ariate fusion/downscaler model, extends the univ ariate downscaler of [ 3 ], and, following [ 3 ], regresses the biv ariate v ector of observed ozone and PM 2 . 5 concen tration on the numerical mo del outputs for both p ollu- tan ts using spatially-v arying co efficien ts [ 19 , 44 ], mo deled as a correlated six-dimensional Gaussian pro cess. The mo del sp ecification is general and flexible, while offering feasible computation for a fully mo del-based fusion across space and time. W e explore sev eral differen t special cases, corresp ond- ing to different assumptions on the correlation structure of the tw o p ollu- tan ts. An imp ortan t feature of our mo del is that it allo ws us to handle not only the spatial misalignment b et w een monitoring data and model output, but also it allo ws us to accommodate the spatial and temp oral misalign- men t b et ween the ozone and PM 2 . 5 monitoring data. Finally , exploiting the correlation, not only b et w een the mo del output and the corresp onding mon- itoring data for a p ollutan t, but also betw een p ollutan ts, we can join tly predict ozone and PM 2 . 5 at any site in the spatial domain and pro vide a measure of the uncertainties asso ciated with suc h predictions. Additionally , at sites where one p ollutan t is measured but the other is not, we can predict the lev el of the latter. The pap er is organized as follo ws. In Section 2, we present the moni- toring data and the numerical model output, in Section 3 we describe the biv ariate do wnscaler mo del, b y first in tro ducing the downscaler in the uni- v ariate setting, and then extending the mo del to the biv ariate case, first in a purely spatial setting and then in a spatio-temp oral setting. Details on ho w to handle the spatio-temp oral misalignment in the monitoring data are also discussed in Section 3. Section 4 presen ts details on the mo del fitting, while Section 5 displa ys results of our analysis. Finally , w e conclude in Section 6 with a discussion and ev aluation of our metho d and with indications for future w ork. 2. Data. Particulate matter (PM 2 . 5 ) and ozone are tw o of the “criteria p ollutan ts” that the Environmen tal Protection Agency (EP A) is required to monitor b y the Clean Air Act. The first is a mixture of solid and liquid particles, emitted in the atmosphere either directly from a source or as result of complicated reactions of chemicals, while the second is a gas made of three atoms of oxygen that is created b y a chemical reaction betw een o xides of nitrogen and v olatile organic comp ounds in the presence of sunlight. The EP A tracks b oth criteria pollutants using both measurements taken at monitoring sites and estimates of air p ollutan ts concen tration pro duced by the Models-3/Communit y Mesoscale Air Qualit y (CMA Q) mo del [ 6 ]( epa.gov/asmdnerl/CMAQ ). imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 5 The latter is a deterministic numerical mo del that predicts concentrations for v arious p ollutan ts b y in tegrating three ma jor components: an atmospheric comp onen t accounting for the atmosphere and its states and motions, an emission comp onen t accoun ting for the emissions injected in the atmosphere, and a c hemical comp onen t accoun ting for the reactions b etw een the different gas present in the atmosphere. By sim ulating v arious chemical and ph ysi- cal pro cesses, suc h as horizon tal and v ertical advection, emission injection, dep osition, plume chemistry effects, etc., CMAQ pro duces estimates of p ol- lutan ts concen tration at pre-determined spatial scales, e.g., 36 km, 12 km, and, most recen tly , at 4 km grid cells. W e consider daily estimates and measurements of ozone and PM 2 . 5 con- cen tration for the region in the South Eastern part of the United States sho wn in Figure 1 during the p erio d June 1 - September 30, 2002, the sum- mer season when solar radiation and temperature are high. In turn this en- courages not only high concentrations of ozone but also of particulate SO 4 , the dominan t component of PM 2 . 5 in the eastern U.S as well as particulate ammonium nitrate. Relativ ely strong association b et ween ozone and PM 2 . 5 concen trations is exp ected and, in fact, observed. F ollowing the air qualit y standards (NAAQS; http://www.epa.go v/air/criteria.h tml) set by EP A for ozone and PM 2 . 5 , daily ozone concen tration is measured as the daily maxi- m um 8-hour a v erage concentration, while daily PM 2 . 5 concen tration is giv en b y the 24-hour a v erage concen tration of particulate matter. The monitoring data used in our case study comes from 226 sites sparsely lo cated in the region, and it has b een obtained from the EP A Air Qual- it y System (AQS) rep ository database. Of these 226 sites, not all measure b oth p ollutan ts: 71 rep ort only ozone measuremen ts, 50 rep ort only PM 2 . 5 , and 105 measure b oth, pinpointing the spatial misalignmen t in the monitor- ing data. In order to v alidate the mo del with out-of-sample predictions, w e randomly split the sites in to t w o sets: a training set used to fit the mo del, and a v alidation set used to assess the p erformance of the model. In the training set, 52 sites measured only ozone concentration, 39 rep orted only PM 2 . 5 concen tration, and 70 rep orted both. In total, the dataset used to fit the mo del comprised 14,630 daily measuremen ts of ozone concen tration and 4,790 measurements of PM 2 . 5 concen tration. The smaller n um b er of obser- v ations for particulate matter is due to the sampling sc heme for PM 2 . 5 : of the 109 PM 2 . 5 monitoring sites, only 11 (10.1%) measure concentration of particulate matter every da y , while 90 (82.5%) rep ort measurements ev ery three da ys, and 8 (7.3%) measure PM 2 . 5 ev ery six da ys. The difference in sampling frequency b et w een ozone and PM 2 . 5 suggests that, in addition to spatial misalignment, there is also temp oral misalign- imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 6 V. J. BERR OCAL ET AL. men t in the monitoring data: ev en sites that measure both p ollutan ts do not sample the t wo p ollutan ts with the same temp oral frequency . Ho w ev er, when concen tration measuremen ts for b oth p ollutan t are av ailable, they display considerable correlation, equal to 0.69. Spatial and temp oral misalignment is also presen t among the v alidation sites: 54 measured only ozone, 46 only PM 2 . 5 sites, while 35 rep orted concentration for b oth p ollutan ts. In total, the v alidation dataset contained 6,530 daily observ ations of ozone concen- tration and 2,559 daily observ ations of concentration of particulate matter. The lo cation of the sites used to fit the and v alidate the model is sho wn in Figure 1 . An exploratory analysis of the observ ed daily concentration data for ozone and PM 2 . 5 rev ealed the need for a transformation. W e modeled the daily ozone concen tration data on the square ro ot scale to stabilize the v ariance, while w e log-transformed the PM 2 . 5 concen tration data to ac hieve normal- it y . Both transformations hav e b een previously used; see, e.g., [ 41 ], [ 40 ], [ 7 ], [ 25 ]. A normal Q-Q plot of the square ro ot of the observ ed daily concen tra- tion of ozone at sites used to fit the mo del is sho wn in Figure 2 (a), while Figure 2 (b) presen ts a normal Q-Q plot of the logarithm of daily a v erage concen tration of PM 2 . 5 at the 109 monitoring sites in the test set measur- ing particulate matter. Both plots seems to suggest a sligh t deviation from normalit y , how ev er the deviation is minimal and we are comfortable using a normal model, in agreement with the literature. Nev erthless, in Section 6 w e address distributional issues more at lengths and suggest an extension of our downscaling approach for non-Gaussian v ariables. The ov erall mean and standard deviation for the square ro ot of ozone concen tration and for the log of PM 2 . 5 concen tration at the test sites w ere, respectively , 7.41 and 1.31 √ ppb (ppb:parts per billions), and 2.72 and 0.56 log( µ g/m 3 ). Estimates of the daily av erage and standard deviation for ozone and PM 2 . 5 are presen ted in Figure 3 . Both panels in Figure 3 indicate daily v ariabilit y in the concen tration lev el for the p ollutan ts, with Figure 3 (a) rev ealing the seasonalit y of ozone during the study p eriod (June 1 - September 20, 2002). In terms of v ariabilit y , b oth ozone and PM 2 . 5 presen t the largest standard deviation on June 25, 2002. F or this day , plots of the observed ozone and PM 2 . 5 concen tration are shown, on the original scales, resp ectiv ely , in Fig- ure 5 (a) and Figure 6 (a). (These plots are developed in conjunction with the data analysis of Section 5.) The n umerical mo del output data consists of estimates of ozone and PM 2 . 5 concen tration lev el generated b y the CMAQ numerical mo del run at 12 km grid cell resolution. An example of the t yp e of concen tration surfaces yielded b y CMA Q can b e observ ed, respectively , in Figure 5 (b) and Figure 6 (b). imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 7 CMA Q produces estimates of a pollutant concentration in terms of a v erage o v er a grid cell (in this case, ov er 10,504 12-km grid cells), while observ a- tions are taken at points. Therefore, there is a spatial misalignmen t b et w een observ ational and numerical model data. W e solv e this difference in spatial resolution b et w een the monitoring and the numerical model data by asso- ciating to eac h site s in the spatial domain S , the CMA Q grid cell B in whic h s lies. This allo ws a comparison b etw een the mo del output and the monitoring data. T o reveal the need for calibration in the numerical model output, w e hav e pro duced pairwise scatterplots of the square ro ot of the observ ed ozone concen tration v ersus the square ro ot of the CMA Q predicted concentration lev el of ozone during the perio d June 1 - Septem b er 30, 2002, of the square ro ot of the observed ozone concen tration versus the log of the CMAQ pre- dicted PM 2 . 5 concen tration, and analogous plots for the log of the observ ed concen tration of particulate matter. W e see that the predictions by the nu- merical mo del are biased and need to b e calibrated; how ev er, they do contain useful information to improv e prediction for b oth pollutants. Additionally , the p ositiv e and substan tial correlation b et ween the CMA Q mo del output for PM 2 . 5 and the observed ozone concentration (r=0.62), and, similarly , b et w een the CMAQ model output for ozone and the observ ed PM 2 . 5 con- cen tration (r=0.69), indicate that the CMA Q output for PM 2 . 5 migh t b e useful to predict ozone and conv ersely . The t ype of calibration implied b y the pairwise scatterplots men tioned ab o v e is constan t across the en tire spatial domain S in consideration. Em- pirical evidence suggests instead that the error and the bias of the numerical mo del output migh t not b e constan t in space, rather they v ary from site to site. Figure 4 presen ts spatial maps of the estimates of the intercept and of the regression co efficien ts of CMA Q ozone and CMAQ PM 2 . 5 at each of the sites used to fit our biv ariate do wnscaler mo del. These estimates hav e b een obtained by regressing at each site s , resp ectiv ely , the normalized log of the observ ed PM 2 . 5 concen tration at s across time on the corresponding normalized CMAQ mo del output for ozone and PM 2 . 5 , b oth taken on the appropriate scale (square ro ot for ozone, log for PM 2 . 5 ). As the plots indi- cate, there is spatial v ariability in the estimates of these coefficients with differen tial v ariabilit y across the estimates. In particular, we see less spatial v ariability in the estimates obtained for the observ ed log PM 2 . 5 concen tra- tion. Also, there is a difference b et ween the t w o p ollutan ts in the significance of the estimates of the co efficien ts: while for ozone, the estimates of the coef- ficien t of CMAQ ozone are all significan t, and ab out 60% of the estimates of the intercept are significant, in the regression for PM 2 . 5 , only 21% of the in- imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 8 V. J. BERR OCAL ET AL. tercepts are significantly differen t from zero, and ab out 88% of the estimates of the coefficient of CMAQ PM 2 . 5 are significan tly differen t from zero. 3. Mo deling. Here w e present the model underlying our downscaling tec hnique. W e first briefly review the univ ariate do wnscaler, presented in [ 3 ], then we extend it to the biv ariate setting. In doing so, we first present the general biv ariate do wnscaler in a purely spatial setting, then we examine extension to accommodate data collected o ver time. 3.1. Univariate downsc aler. Let Y ( s ) b e the observed data at s and let x ( B ) b e the n umerical model output at grid cell B . Then, x ( B ) is inter- preted as the av erage lev el of the v ariable under consideration o v er B . The spatial downscaler addresses the difference in spatial resolution b et ween the observ ational data and the numerical mo del output, i.e. the “change of sup- p ort” problem [ 2 , 9 ], b y associating to eac h p oin t s the grid cell B in whic h s lies. Then, it relates the observ ational data and the numerical model output as follo ws: (1) Y ( s ) = ˜ β 0 ( s ) + ˜ β 1 ( s ) x ( B ) +  ( s )  ( s ) iid ∼ N (0 , τ 2 ) where ˜ β 0 ( s ) = β 0 + β 0 ( s ) ˜ β 1 ( s ) = β 1 + β 1 ( s ) (2) In ( 2 ), β 0 and β 1 denote, respectively , the o v erall additiv e and m ultiplicativ e bias of the numerical mo del, and β 0 ( s ) and β 1 ( s ) are lo cal adjustmen ts to the ov erall bias terms. Finally ,  ( s ) is a white noise pro cess with nugget v ariance τ 2 . Eviden tly , a constan t-mean downscaler is a sp ecial case of ( 2 ) obtained when the lo cal adjustmen ts, β 0 ( s ) and β 1 ( s ), are set to zero. An ticipating association b et w een intercept and slop e, at the second hier- arc hical lev el of the mo del, w e set the t w o spatially-v arying co efficien ts β 0 ( s ) and β 1 ( s ) to b e a biv ariate correlated mean-zero Gaussian pro cesses. Using the metho d of coregionalization [ 19 , 44 , 51 ], we assume that there exist t wo laten t mean-zero, unit-v ariance indep endent Gaussian pro cesses w 0 ( s ) and w 1 ( s ) such that Co v( w j ( s ) , w j ( s 0 )) = exp( − φ j k s − s 0 k ), where, for j = 0 , 1, φ j is the spatial deca y parameter of the Gaussian process w j ( s ), k s − s 0 k is imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 9 the great circle distance b et ween s and s 0 1 , and (3) β 0 ( s ) β 1 ( s ) ! = A w 0 ( s ) w 1 ( s ) ! The coregionalization matrix A in ( 3 ) determines the correlation b et w een the t w o spatially-v arying coefficients β 0 ( s ) and β 1 ( s ) and can b e assumed, without loss of generalit y , to b e low er-triangular. Note that T = AA 0 is the constant lo cal cov ariance matrix. The sp ecification of the mo del is then completed with the priors for: the ov erall bias terms, β 0 and β 1 , the nugget v ariance τ 2 , the three non-n ull elemen ts of the coregionalization matrix A , and the decay parameters φ 0 and φ 1 . Other choices for obtaining dep enden t Gaussian pro cesses could b e entertained, including moving a verages [ 26 ], con v olution of co v ariance functions [ 34 ], and particular parametric c hoices. In the last case, some p ossibilities include the cross-co v ariance functions in tro duced by [ 1 ] and [ 20 ]. 3.2. The bivariate downsc aler: static version. T urning to the biv ariate case, we describ e the mo del in terms of square ro ot ozone and log PM 2 . 5 concen trations, recognizing that the model could b e applied to an y suitable pair of point referenced v ariables. Let Y 1 ( s ) and Y 2 ( s ) denote, respectively , the square ro ot and the logarithm of the observ ed ozone and PM 2 . 5 con- cen tration at s , and let x 1 ( B ) and x 2 ( B ) indicate, respectively , the square ro ot and the logarithm of the n umerical mo del output for ozone and PM 2 . 5 at grid cell B . Then, Y ( s ) = ( Y 1 ( s ) , Y 2 ( s )) is the biv ariate random vec- tor with the observ ed concentrations for the tw o p ollutants at s , while x ( B ) = ( x 1 ( B ) , x 2 ( B )) is the biv ariate v ector with the t wo n umerical mo del outputs. Building upon ( 1 ), our biv ariate do wnscaler model relates the mon- itoring data at s and the n umerical mo del outputs at grid cell B , with s lying in B , as follo ws: Y 1 ( s ) = ˜ β 10 ( s ) + ˜ β 11 ( s ) x 1 ( B ) + ˜ β 12 ( s ) x 2 ( B ) +  1 ( s ) Y 2 ( s ) = ˜ β 20 ( s ) + ˜ β 21 ( s ) x 1 ( B ) + ˜ β 22 ( s ) x 2 ( B ) +  2 ( s ) (4) where  1 ( s ) and  2 ( s ) are tw o independent white noise processes with n ugget v ariances, resp ectiv ely , τ 2 1 and τ 2 2 . 1 When considering large spatial domains, it is preferable to use three-dimensional Eu- clidean distance as an argumen t for the exponential co v ariance function rather than great- circle distance, which migh t yield a co v ariance matrix which is not positive definite. F or the domain w e study , the tw o metrics almost coincide and yield v ery similar results. imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 10 V. J. BERR OCAL ET AL. As in the univ ariate case, we adopt a random effect notation, decomp osing eac h of the β ij ( s ), i = 1 , 2, j = 0 , 1 , 2, in to the sum of an o v erall term and a lo cal adjustment. That is, (5) ˜ β ij ( s ) = β ij + β ij ( s ) . The six-dimensional pro cess ( β ij ( s )) i =1 , 2; j =0 , 1 , 2 is in turn modeled as a six- dimensional mean-zero correlated Gaussian pro cess, again using the method of coregionalization. Therefore, w e express eac h of the spatially-v arying co- efficien ts β ij ( s ) as a linear combination of mean-zero unit-v ariance laten t indep enden t Gaussian pro cesses, each equipp ed with an exp onen tial co v ari- ance structure. More sp ecifically , assuming, without loss of generalit y , that the coregionalization matrix A is low er triangular, we hav e that: (6)          β 10 ( s ) β 11 ( s ) β 12 ( s ) β 20 ( s ) β 21 ( s ) β 22 ( s )          = A          w 1 ( s ) w 2 ( s ) w 3 ( s ) w 4 ( s ) w 5 ( s ) w 6 ( s )          where, for eac h k = 1 , . . . , 6, Cov( w k ( s ) , w k ( s 0 )) = exp( − φ k k s − s 0 k ), with φ k deca y parameter for w k ( s ). As in the univ ariate case, the coregionalization matrix A determines, through T = AA 0 , the correlation b et w een the six spatially-v arying co efficients β ij ( s ). Additionally , as a consequence of ( 4 ), A induces a correlation structure on the biv ariate random vector Y ( s ), i.e. (7) Σ Y ( s ) ,Y ( s 0 ) =  I 2 ⊗ (1 x 1 ( B ) x 2 ( B )) 0  · AΣ w A 0 ·  I 2 ⊗ (1 x 1 ( B 0 ) x 2 ( B 0 ))  where B and B 0 are grid cells con taining, resp ectively , s and s 0 , I 2 is the iden tit y matrix of order 2, and Σ w is a 6 × 6 diagonal matrix with i -th elemen t exp( − φ i k s − s 0 k ), i = 1 , . . . , 6. Our biv ariate do wnscaler mo del, though simple, is general and flexible. The sp ecification of the bias of the numerical mo del by means of spatially- v arying co efficien ts recognizes that calibration under the n umerical mo del should b e site-sp ecific. Again, as in the univ ariate case, a constant-mean bi- v ariate downscaler can b e obtained as a special case. Moreo v er, p erm utation of the en tries in the β ( s ) vector do es not affect the prior. The joint mo del still presents a 6 × 6 lo cal cov ariance matrix mo deled through its Cholesky decomp osition. Simplified v ersions of the general biv ariate do wnscaler can b e obtained b y imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 11 setting entries in the coregionalization matrix A equal to zero. 2 F or example, if the set of non-n ull entries in the A matrix is giv en b y { A 11 , A 21 , A 22 , A 44 , A 64 , A 66 } and w e assume that β 12 ≡ 0 and β 21 ≡ 0, then our biv ariate do wnscaler reduces to t w o indep enden t univ ariate downscalers on ozone and PM 2 . 5 , re- sp ectiv ely . This form of A implies that the spatially-v arying co efficien ts are correlated within pollutants but not across p ollutan ts. A simpler v ersion of the biv ariate downscaler that requires only three non-n ull entries in the coregionalization matrix A but still induces a cor- relation among the tw o comp onen ts Y 1 ( s ) and Y 2 ( s ) of Y ( s ) is the mo del that assumes that only the lo cal intercept adjustmen ts, β 10 ( s ) and β 20 ( s ), are non-null and correlated. In this case, the coregionalization matrix A has only three non-n ull en tries, { A 11 , A 41 , A 44 } and the co v ariance b et ween the square ro ot of the observ ed ozone concen tration at s and the logarithm of the observ ed PM 2 . 5 concen tration at s 0 reduces simply to: (8) Co v( Y 1 ( s ) , Y 2 ( s 0 )) = A 11 · A 41 exp( − φ 1 k s − s 0 k ) A simple extension to this mo del can b e ac hieved by maintaining the same correlation structure for the biv ariate random v ector Y ( s ), that is, b y assuming again that A 41 is the only non-null off-diagonal elemen t of the core- gionalization matrix A , but by p ostulating that the lo cal adjustmen ts to the co efficien ts of the tw o numerical mo del outputs, β 11 ( s ), β 12 ( s ), β 21 ( s ), and β 22 ( s ), are non-null indep enden t Gaussian pro cesses with v ariances equal, resp ectiv ely , to A 2 22 , A 2 33 , A 2 55 and A 2 66 . Then, the coregionalization matrix A corresponding to this mo del has seven non-null entries: A 11 , A 22 , A 33 , A 41 , A 44 , A 55 , and A 66 . An extension of this last mo del that still inv olv es fewer parameters than the general form ulation for the biv ariate downscaler (13 v ersus 21) is the mo del which assumes that only lo cal adjustments of the form β 1 j ( s ) and β 2 j ( s ), for j = 0 , 1 , 2, are correlated across p ollutan ts, while all the β ij ( s )’s sp ecific to the same p ollutan t are correlated. Such a mo del sp ecification corresp onds to the following coregionalization matrix: (9) A =          A 11 0 0 0 0 0 A 21 A 22 0 0 0 0 A 31 0 A 33 0 0 0 A 41 0 0 A 44 0 0 0 A 52 0 A 54 A 55 0 0 0 A 63 A 64 0 A 66          2 In principle, we could include all 21 A ij ’s in the mo del and let the data suggest which are significant. Instead, w e ha ve c hosen to do model comparison betw een several mo dels, eac h having plausible interpretations. imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 12 V. J. BERR OCAL ET AL. Other simplifications of the biv ariate downscaler general mo del can be en- visioned. The sp ecification of our Ba y esian hierarc hical biv ariate downscaler mo del is completed with priors on the parameters. W e adopt an inv erse gamma with large v ariance for the nugget v ariances, τ 2 1 and τ 2 2 ; indep enden t v ague normals for the o verall regression co efficien ts β ij , i = 1 , 2, j = 0 , 1 , 2; log- normals with v ague standard deviations for the diagonal entries of the core- gionalization matrix A ; and v ague normals for the off-diagonal elemen ts of A . Details on handling the decay parameters φ k , k = 1 , . . . , 6, are pro vided in Section 4.1 . 3.3. Sp ac e-time bivariate downsc aler. W e no w extend our general biv ari- ate downscaler to accommo date data that has b een collected ov er time. Let t = 1 , . . . , T denote the times at whic h we hav e observ ations and n umerical mo del outputs, and let Y 1 ( s , t ) and Y 2 ( s , t ) denote, respectively , the square ro ot and the logarithm of the observ ed ozone and PM 2 . 5 concen tration at site s on da y t . In an analogous wa y , let x 1 ( B , t ) and x 2 ( B , t ) indicate, re- sp ectiv ely , the square root and the logarithm of the CMAQ output for ozone and PM 2 . 5 on da y t o ver the 12-km grid cell B . Then, if s lies in grid cell B , the mo del for our biv ariate downscaler b ecomes Y 1 ( s , t ) = ˜ β 10 ( s , t ) + ˜ β 11 ( s , t ) x 1 ( B , t ) + ˜ β 12 ( s , t ) x 2 ( B , t ) +  1 ( s , t ) Y 2 ( s , t ) = ˜ β 20 ( s , t ) + ˜ β 21 ( s , t ) x 1 ( B , t ) + ˜ β 22 ( s , t ) x 2 ( B , t ) +  2 ( s , t ) (10) where  1 ( s , t ) and  2 ( s , t ) are tw o white noise processes that follo w indep en- den tly a N (0 , τ 2 1 ) and a N (0 , τ 2 2 ) distribution. As in the spatial setting, for each i = 1 , 2 and j = 0 , 1 , 2, w e write: (11) ˜ β ij ( s , t ) = β ij,t + β ij ( s , t ) where, for eac h t = 1 , . . . , T , β ij ( s , t ) are correlated Gaussian pro cesses. F ollowing [ 3 ], w e can model the temp oral dep endence in the comp onen ts of ( 11 ) in several w ays. F or example, for the β ij,t ’s, we can assume that they are either indep enden t across time, i.e. β ij,t ind ∼ N ( µ ij , σ 2 ij ), or, alternatively , that they ev olv e dynamically in time [ 52 ], that is, (12) β ij,t = ρ ij β ij,t − 1 + η ij,t , η ij,t ind ∼ N (0 , ξ 2 ij ) , i = 1 , 2; j = 0 , 1 , 2 and β ij, 0 ∼ N ( µ ij, 0 , σ 2 ij, 0 ). Similar time-dependence can be imp osed on the Gaussian processes β ij ( s , t ). W e can assume that the correlated Gaussian pro cesses β ij ( s , t ) are of the imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 13 form (13)          β 10 ( s , t ) β 11 ( s , t ) β 12 ( s , t ) β 20 ( s , t ) β 21 ( s , t ) β 22 ( s , t )          = A          w 1 ( s , t ) w 2 ( s , t ) w 3 ( s , t ) w 4 ( s , t ) w 5 ( s , t ) w 6 ( s , t )          with A still a coregionalization matrix and with the underlying Gaussian pro cesses, w k ( s , t ), nested within time. Then, at each t = 1 , . . . , T , the w k ( s , t ) are independent replicates of mean-zero unit-v ariance Gaussian pro- cesses with exp onen tial correlation and deca y parameters φ k,t that can b e either tak en to be constant in time or indep enden t across time. In the second case, w e can mo del the lo cal adjustments, β ij ( s , t ), to evolv e dynamically in time. As in [ 17 ], for eac h i = 1 , 2 and j = 0 , 1 , 2, we assume that (14) β ij ( s , t ) = γ ij β ij ( s , t − 1) + ν ij ( s , t ) where the innov ations ν ij ( s , t ) are correlated Gaussian processes. In other w ords, in this second case, the coregionalization ( 13 ) is sp ecified on the ν ij ( s , t ), rather than on the β ij ( s , t ), but it still employs the underlying Gaus- sian pro cesses w k ( s , t ), which are defined as abov e. This dynamic mo del for the β ij ( s , t ) is completed b y sp ecifying as initial conditions that β ij ( s , 0) = 0 for i = 1 , 2 and j = 0 , 1 , 2. In b oth cases w e can b e more general and assume that the coregionalization matrix in ( 13 ) is indexed by time, with entries indep enden t across time. F urthermore, if t is on a con tinuous domain, then the w k ( s , t ) can be indep enden t space time processes. The tw o different time dep endence structures for the o v erall regression co efficien ts, β ij,t , and for their lo cal adjustments, β ij ( s , t ), can b e combined together in four wa ys, yielding four mo dels that corresp ond to different as- sumptions on the wa y the o v erall and local performance and calibration of the n umerical mo del changes ov er time. 3.4. The gener al multivariate downsc aler. Extension to a m ultiv ariate do wnscaler given observ ational data and numerical model outputs for p ran- dom v ariables is evident. As in the biv ariate case, w e start with the static form ulation of the model and then w e extend it to the spatio-temp oral set- ting. Let Y i ( s ), i = 1 , . . . , p b e the observ ed data for the i -th v ariable at a site s in the spatial domain S and let x i ( B ), i = 1 , . . . , p , be the n umerical mo del imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 14 V. J. BERR OCAL ET AL. output for the i -th v ariable ov er grid cell B . Then, again w e asso ciate to eac h site s the grid cell B in whic h s lies and we relate the observ ational data and the numerical model output as follo ws: (15) Y i ( s ) = ˜ β i 0 ( s ) + p X j =1 ˜ β ij ( s ) x j ( B ) +  i ( s )  i ( s ) iid ∼ N (0 , τ 2 i ) for i = 1 , . . . , p . W e decomp ose eac h of the p · ( p + 1) terms ˜ β ij ( s ), i = 1 , . . . , p , j = 0 , 1 , . . . , p in the sum of an o v erall term and a lo cal adjustment: ˜ β ij ( s ) = β ij + β ij ( s ), and w e mo del the β ij ( s )’s as correlated mean-zero Gaussian pro cesses with exp onen tial cov ariance structures using again the metho d of coregionalization. In the general p -dimensional multiv ariate do wnscaler, the coregionalization matrix A is a p ( p + 1) × p ( p + 1) matrix. Sp ecifications of ( 15 ) similar to the ab ov e can b e considered by simply setting to zero en tries of the A matrix, thus inducing simplifications on b oth the correlation structure of the m ultiv ariate random v ector Y = { ( Y i ( s )) i =1 ,...,p : s ∈ S } and on the co v ariance structure of the single Y i ( s ), i = 1 . . . , p . T o complete the specification of the mo del, we place standard priors on the mo del parameters: thus, w e sp ecify p ( p + 1) v ague normals for the o verall regression co efficien ts β ij ; p in v erse gammas with large v ariances for the n ugget v ariances τ 2 i , i = 1 . . . , p ; log-normals with large v ariances for the diagonal entries A kk , k = 1 , . . . , p ( p + 1), of the coregionalization matrix A ; and v ague normals for the off-diagonal entries of A . If the data on the p random v ariables has b een collected not only ov er space, but also ov er time, then the m ultiv ariate general downscaler would b e extended to a spatio-temp oral setting in an ob vious wa y . Let Y i ( s , t ), i = 1 , . . . , p b e the observ ed data for the i -th random v ariable at site s at time t , t = 1 , . . . , T , and let x i ( B , t ) be the mo del output for the i -th v ariable o v er grid cell B at time t . Then, for i = 1 , . . . , p , we p ostulate the following relationship b et w een observ ational data and numerical model output: (16) Y i ( s , t ) = ˜ β i 0 ( s , t ) + p X j =1 ˜ β ij ( s , t ) x j ( B , t ) +  i ( s , t )  i ( s , t ) iid ∼ N (0 , τ 2 i ) where s lies in grid cell B . F or each i = 1 , . . . , p and j = 0 , 1 , . . . , p , we write: ˜ β ij ( s , t ) = β ij,t + β ij ( s , t ), and we mo del the temp oral structure in the β ij,t and in the β ij ( s , t ), follo wing what we hav e presented in Section 3.3 for the biv ariate do wnscaler. Th us, we assume that the β ij,t ’s are either nested within time or dynamic in time, and, similarly , we p ostulate that the β ij ( s , t ) are either indep enden t replicates o v er time or are dynamic in time. imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 15 4. Mo del fitting. 4.1. Priors. W e ha v e already men tioned briefly in Sections 3.2 and 3.4 the form of the prior distributions for the parameters of the biv ariate and m ultiv ariate do wnscaler both in their static and spatio-temporal form ula- tions. Here, we turn to the spatial deca y parameters φ k of the latent Gaussian pro cesses w k ( s ) used in the coregionalization. Discrete uniform priors facili- tate mo del fitting. Previous exp erience with such priors has alwa ys resulted in p osterior distributions placing highest probabilit y on a v alue which is, essen tially , the Restricted Maximum Lik eliho od (REML; [ 24 , 38 ]) estimate. Hence, in what follo ws, w e prop ose to keep the spatial decay parameters fixed and equal to v alues that are determined by a sensitivity analysis. More sp ecifically , consider the simplified biv ariate downscaler that uses only the tw o laten t Gaussian pro cesses w 1 ( s ) and w 4 ( s ) and is obtained when the coregionalization matrix A has only three non-null entries, A 11 , A 41 , and A 44 . The deca y parameters φ 1 and φ 4 determine the w a y the spatial correlation decays with distance in w 1 ( s ) and w 4 ( s ), resp ectiv ely , and, b y consequence, in Y 1 ( s ) and Y 2 ( s ). In fact, the cov ariance structure of Y 2 ( s ) induced by such simplified version of the biv ariate do wnscaler is given b y the sum of tw o exp onen tial cov ariance functions with deca y parameters, resp ectiv ely , equal to φ 1 and φ 4 , and a diagonal matrix with en tries all equal to the nugget v ariance τ 2 2 . T o obtain a rough estimate of the magnitude of the deca y parameters φ 1 and φ 4 , we hav e pro ceeded as follo ws. W e ha ve considered the moni- toring data for the square ro ot of ozone and the logarithm of PM 2 . 5 for a giv en day t . W e hav e modeled the monitoring data for b oth v ariables as Gaussian pro cesses with an unknown mean, resp ectively , µ Y 1 ,t and µ Y 2 ,t , and with an exponential co v ariance function plus nugget effect, with decay parameters φ Y 1 ,t and φ Y 2 ,t , and n ugget v ariance, τ 2 Y 1 ,t and τ 2 Y 2 ,t . W e ha v e esti- mated the parameters via Restricted Maximum Likelihoo d (REML; [ 24 , 38 ] as implemented in the geoR pac k age of R . W e ha v e rep eated this op eration for eac h da y in the p eriod June 1– Septem b er 30, 2002 and ha v e obtained daily estimates of φ Y 1 ,t and φ Y 2 ,t . Since histograms of those daily estimates exhibited long right tails, we ha v e summarized those distributions by tak- ing the median of the daily estimates, whic h yielded, resp ectiv ely , 0 . 0016 and 0 . 00125, corresp onding to ranges of, respectively , 1875 and 2400 km. W e ha v e then p erformed a sensitivity analysis to determine ho w the pre- dictiv e performance of a downscaler mo del is affected b y the magnitude of the decay parameter. Therefore, keeping the decay parameters φ 1 and φ 4 fixed and equal to the estimated medians of φ Y 1 ,t and φ Y 2 ,t , w e hav e fit the imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 16 V. J. BERR OCAL ET AL. spatio-temp oral v ersion of the biv ariate downscaler that models the ov erall terms β ij,t , i = 1 , 2, j = 0 , 1 , 2 and the lo cal adjustmen ts to the o v erall in tercepts, β i 0 ( s , t ), i = 1 , 2 as indep endent in time (see Section 7.2 in the Supplemen tal material online). F or each day , w e ha v e predicted levels of ozone and PM 2 . 5 concen tration at the v alidation sites s 0 b y sampling from the p osterior predictiv e distributions, f ( Y ( s 0 ) | { Y ( s ) } , { x ( B ) } ) and subse- quen tly backtransforming the prediction to the original scale. W e hav e then compared the predictions with the observ ations by computing the predictive mean square error (PMSE), predictiv e mean absolute error (PMAE), empir- ical co verage of the 95% predictive interv al and width of the 95% predictive in terv al. More sp ecifically , denoting with Y ? ( s , t ) the bac ktransformed data (recall w e used the square ro ot and the logarithm transform on ozone and PM 2 . 5 , resp ectiv ely), we hav e computed the PMSE as follows: (17) PMSE = 1 n v T X t =1 V X r =1  c Y ? ( s r , t ) − Y ? ( s r , t )  2 · I ( Y ? ( s r , t )) where n v denotes the total num b er of observ ations in the v alidation dataset (6530 for ozone or 2559 for PM 2 . 5 as noted in Section 2 ), s r is the r -th site in the v alidation set con taining a total of V sites, c Y ? ( s r , t ) denotes the p osterior mean of the predictive distribution at site s r at time t on the original scale, and I ( Y ? ( s r , t )) is equal to 1 if Y ? ( s r , t ) is observ ed and 0 otherwise. An analogous definition holds for the PMAE, with c Y ? ( s r , t ) no w referring to the p osterior median of the predictiv e distributions, again on the original scale. W e ha ve generated predictions and computed the summary statistics men tioned abov e four times, each time keeping the deca y parameters fixed and setting them equal, resp ectively , to the estimated medians of φ Y 1 ,t and φ Y 2 ,t , m ultiplied b y 10 and 100, and divided by 10 and 100. T o compare the performance of the biv ariate do wnscaler to that of a uni- v ariate downscaler, we hav e p erformed the same sensitivit y analysis also for the univ ariate do wnscaler. Thus, keeping the deca y parameters for ozone and PM 2 . 5 fixed and equal to the v alues considered for the biv ariate do wn- scaler, we ha ve fit tw o spatio-temporal univ ariate indep endent do wnscalers for ozone and particulate matter, where the only non-null local adjustment term w as the adjustment to the o v erall intercept, and where b oth the o ver- all regression terms and the spatially v arying co efficien ts w ere nested within time. Then, as in the biv ariate do wnscaler case, we hav e predicted ozone and PM 2 . 5 concen tration at the v alidation sites and we hav e assessed the p er- formance of the predictions using the same summary statistics mentioned ab o v e. W e hav e carried out a similar sensitivit y analysis also for ordinary kriging imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 17 [ 8 , 9 ]. Therefore, for b oth ozone and PM 2 . 5 , for eac h da y , using the same range of v alues for the deca y parameters used in the experiments with the univ ariate and biv ariate do wnscalers, w e hav e kriged the observ ed data, on the transformed scale, to the v alidation sites, and we hav e subsequen tly bac k- transformed the predictions. Finally , w e ha v e compared the predictions with the observ ations using the same summary statistics that w e ha v e emplo yed for b oth downscalers. T able 1 and T able 2 rep ort, resp ectively , the Predictive Mean Square Er- ror (PMSE), Predictiv e Mean Absolute Error (PMAE), the empirical cov- erage and the width of the nominal 95% predictiv e interv al for ozone and PM 2 . 5 for the tw o space-time do wnscalers and for ordinary kriging. In both tables, the results corresp onding to predictions obtained when the mo dels ha v e b een fitted using v alues of the decay parameters equal to the medians of the daily estimates of φ Y 1 ,t and φ Y 2 ,t are rep orted in the middle column. F rom both tables, it is clear that when the decay parameter is one and, esp ecially , tw o orders of magnitude larger, the qualit y of the predictions ob- tained using b oth downscalers deteriorates noticeably in terms of PMSE and PMAE. Also, the predictive interv als are wider in this case, an indication that there is increased v ariability in the predictions. Missp ecifying the mag- nitude of the deca y parameters on the small side affects the qualit y of the predictions less noticeably . So, altogether, the sensitivit y analysis has shown that b oth the univ ariate and the biv ariate do wnscaler and ordinary kriging yield predictions that ha v e o v erall best v alidation statistics when the decay parameters are equal to the median of the daily estimates of φ Y 1 ,t and φ Y 2 ,t and so, in the sequel, we fix them at these v alues. 4.2. Hand ling misalignment. As noted in Section 2 , not all sites report- ing ozone concen tration measure PM 2 . 5 and vice v ersa. This creates a spatial misalignmen t in the data that we handle through the laten t Gaussian pro- cesses w k ( s ) in v olved in the coregionalization, after ha ving appropriately p erm uted and partitioned the data vector. Here, we illustrate how to han- dle the spatial misalignmen t in the static setting. Extension to the spatio- temp oral case is straightforw ard. Let S t b e the set of n t sites rep orting measuremen ts of ozone and/or PM 2 . 5 concen tration on da y t . W e can decomp ose S t as the union of three disjoin t sets, S Both ,t , S O,t , and S P M,t : (18) S t = S Both ,t ∪ S O,t ∪ S P M,t where the first set includes all sites s ∈ S in which b oth ozone and PM 2 . 5 ha v e b een measured on day t , the second contains all the sites in which only imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 18 V. J. BERR OCAL ET AL. ozone was measured on day t , and the last includes all sites s where only PM 2 . 5 has b een measured on da y t . F ollowing the decomp osition in ( 18 ), w e reorder and partition the random v ectors Y = ( Y 1 ( s ) , Y 2 ( s )) s ∈S t and w k = ( w k ( s )) s ∈S t in the follo wing wa y: Y = ( Y Both , Y O , Y P M ) and w k = ( w k, Both , w k,O , w k,P M ) where Y Both = ( Y 1 ( s ) , Y 2 ( s )) s ∈S Both ,t and analogous definitions hold for Y O , Y P M , w k, Both , w k,O , and w k,P M . Then it is clear that b oth comp onen ts Y 1 , Both = ( Y 1 ( s )) s ∈S Both ,t and Y 2 , Both = ( Y 2 ( s )) s ∈S Both ,t ha v e non-missing v alues and b oth con tribute to the log-lik eliho od ( ?? ), while for s ∈ S O,t and s ∈ S P M,t only one comp onen t, resp ectiv ely , Y 1 ,O and Y 2 ,P M , has non-missing v alues and con tributes to the log-lik eliho od. The latent Gaussian processes w k ( s ) are defined on the entire spatial do- main S . T o obtain samples from the p osterior distribution of w k ( s ), within eac h MCMC iteration, w e dra w samples from the full conditionals b y pro- ceeding in the follo wing wa y . If Θ denotes the collection of all mo del pa- rameters (thus, Θ = ( τ 2 1 , τ 2 2 , { A kl } k,l =1 ,..., 6; l ≤ k , { β ij } i =1 , 2; j =0 , 1 , 2 )), for each k = 1 , . . . , 6, w e: i. sample w k, Both from the full conditional π ( w k, Both | w k,O , w k,P M , Y Both , Θ); ii. sample w k,O from the full conditional π ( w k,O | w k, Both , w k,P M , Y 1 ,O , Θ); iii. sample w k,P M from the full conditional π ( w k,P M | w k, Both , w k,O , Y 2 ,P M , Θ); iv. for eac h s ∈ S \ S t , sample w k ( s ) from the conditional distribution π ( w k ( s ) | w k, Both , w k,O , w k,P M , Θ) By rep eating steps i - iv within each MCMC iteration, we obtain samples from the p osterior distribution of w k ( s ) for s ∈ S , thus solving the problem of spatial misalignmen t betw een ozone and PM 2 . 5 . 4.3. Mo del sele ction. Using results from an analysis on the predictiv e p erformance of the univ ariate spatio-temp oral downscaler for ozone concen- tration [ 3 ], we consider only the spatio-temp oral versions of the biv ariate do wnscaler with time-v arying parameters indep enden t across time. More sp ecifically , w e hav e fitted the spatio-temp oral v ersions of the four biv ariate do wnscalers that w e ha v e presen ted in Section 3.2 , that is: i. the biv ariate do wnscaler equiv alen t to tw o indep endent univ ariate down- scalers, obtained when the coregionalization matrix A has as non-null en tries only { A 11 , A 21 , A 22 , A 44 , A 64 , A 66 } ; ii. the biv ariate downscaler with only { A 11 , A 41 , A 44 } as non-n ull en tries in the coregionalization matrix A ; imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: October 29, 2018 A BIV ARIA TE SP ACE-TIME DOWNSCALER 19 iii. the biv ariate do wnscaler with { A 11 , A 22 , A 33 , A 41 , A 44 , A 55 , A 66 } as non- n ull en tries in the coregionalization matrix A ; and iv. the biv ariate downscaler with coregionalization matrix A giv en by ( 9 ) F or each of these mo dels, w e ha ve examined their out-of-sample predic- tiv e p erformance, by predicting daily concentrations of ozone and PM 2 . 5 at the 65 v alidation sites describ ed in Section 2 and shown in Figure 1 . These predictions are compared with the observed ozone and PM 2 . 5 lev- els at the v alidation sites during the p eriod June 1-Septem b er 30, 2002. Note that, for eac h day , predictions of ozone and PM 2 . 5 at a v alidation site s 0 are obtained by sampling from the p osterior predictive distribu- tion f ( Y ( s 0 , t ) | { Y ( s , t ) } , { x ( B , t ) } ) and then by subsequen tly transforming them bac k to the original scale. T o quan tify the predictive p erformance of each biv ariate downscaler w e ha v e utilized the same v alidation statistics used in our sensitivity analysis presen ted in Section 4.1 . W e also employ t wo prop er scoring rules: (i) the Con tin uous Rank ed Probabilit y Score (CRPS; [ 21 ]) defined as: CRPS( F , y ) = Z ( F ( z ) − 1 { z ≥ y } ) 2 d z where F is the cumulativ e predictive distribution function, y is the observ a- tion that materializes, 1 is the Hea viside function, that is equal to 1 if z is greater than y and 0 otherwise, and (ii) the Interv al Score [ 21 ], defined for a (1 − α ) · 100% predictive in terv al with lo w er b ound l and upp er b ound u as: IS( l ; u ; y ) =  ( u − l ) + 2 α ( l − y ) 1 { y

Original Paper

Loading high-quality paper...

Comments & Academic Discussion

Loading comments...

Leave a Comment