1. NRCSE Modelling non-stationarity in space and time for air quality dataPeter GuttorpUniversity of Washingtonpeter@stat.washington.edu

2. Outline • Lecture 1: Geostatistical tools • Gaussian predictions • Kriging and its neighbours • The need for refinement • Lecture 2: Nonstationary covariance estimation • The deformation approach • Other nonstationary models • Extensions to space-time • Lecture 3: Putting it all together • Estimating trends • Prediction of air quality surfaces • Model assessment

3. Research goals in air quality modeling • Create exposure fields for health effects modeling • Assess deterministic air quality models • Interpret environmental standards • Enhance understanding of complex systems

4. The geostatistical setup • Gaussian process • m(s)=EZ(s) Var Z(s) < ∞ • Z is strictly stationary if • Z is weakly stationary if • Z is isotropic if weakly stationary and

5. The problem • Given observations at n sites Z(s1),...,Z(sn) • estimate • Z(s0) (the process at an unobserved site) • or (a weighted average of the process)

6. A Gaussian formula • If • then

7. Simple kriging • Let X = (Z(s1),...,Z(Sn))T, Y = Z(s0), so • mX=m1n, mY=m, • SXX=[C(si-sj)], SYY=C(0), and SYX=[C(si-s0)]. • Thus • This is the best linear unbiased predictor for known m and C (simple kriging). • Variants: ordinary kriging (unknown m) • universal kriging (m=Ab for some covariate A) • Still optimal for known C. • Prediction error is given by

8. The (semi)variogram • Intrinsic stationarity • Weaker assumption (C(0) need not exist) • Kriging can be expressed in terms of variogram

9. Estimation of covariance functions • Method of moments: square of all pairwise differences, smoothed over lag bins • Problem: Not necessarily a valid variogram

10. Least squares • Minimize • Alternatives: • fourth root transformation • weighting by 1/g2 • generalized least squares

11. Fitted variogram

12. Kriging surface

13. Kriging standard error

14. A better combination

15. Maximum likelihood • Z~Nn(m,S) S = a [r(si-sj;q)] = a V(q) • Maximize • and maximizes the profile likelihood

16. A peculiar ml fit

17. Some more fits

18. All together now...

19. Effect of estimating covariance structure • Standard geostatistical practice is to take the covariance as known. When it is estimated, optimality criteria are no longer valid, and “plug-in” estimates of variability are biased downwards. • (Zimmerman and Cressie, 1992) • A Bayesian prediction analysis takes proper account of all sources of uncertainty (Le and Zidek, 1992)

20. Violation of isotropy

21. General setup • Z(x,t) = m(x,t) + n(x)1/2E(x,t) + e(x,t) • trend + smooth + error • We shall assume that m is known or constant • t = 1,...,T indexes temporal replications • E is L2-continuous, mean 0, variance 1, independent of the error e • C(x,y) = Cor(E(x,t),E(y,t)) • D(x,y) = Var(E(x,t)-E(y,t)) (dispersion)

22. Geometric anisotropy • Recall that if we have an isotropic covariance (circular isocorrelation curves). • If for a linear transformation A, we have geometric anisotropy (elliptical isocorrelation curves). • General nonstationary correlation structures are typically locally geometrically anisotropic.

23. The deformation idea • In the geometric anisotropic case, write • where f(x) = Ax. This suggests using a general nonlinear transformation . Usually d=2 or 3. • G-plane D-space • We do not want f to fold.

24. Implementation • Consider observations at sites x1, ...,xn. Let be the empirical covariance between sites xi and xj. Minimize • where J(f) is a penalty for non-smooth transformations, such as the bending energy

25. SARMAP • An ozone monitoring exercise in California, summer of 1990, collected data on some 130 sites.

26. Transformation • This is for hr. 16 in the afternoon

27. Thin-plate splines Linear part

28. A Bayesian implementation • Likelihood: • Prior: • Linear part: • fix two points in the G-D mapping • put a (proper) prior on the remaining two parameters • Posterior computed using Metropolis-Hastings

29. California ozone

30. Posterior samples

31. Other applications • Point process deformation (Jensen & Nielsen, Bernoulli, 2000) • Deformation of brain images (Worseley et al., 1999)

32. Isotropic covariances on the sphere Isotropic covariances on a sphere are of the form where p and q are directions, gpq the angle between them, and Pi the Legendre polynomials. Example: ai=(2i+1)ri

33. A class of global transformations • Iteration between simple parametric deformation of latitude (with parameters changing with longitude) and similar deformations of longitude (changing smoothly with latitude). • (Das, 2000)

34. Three iterations

35. Global temperature • Global Historical Climatology Network 7280 stations with at least 10 years of data. Subset with 839 stations with data 1950-1991 selected.

36. Isotropic correlations

37. Deformation

38. Assessing uncertainty

39. Gaussian moving averages • Higdon (1998), Swall (2000): • Let x be a Brownian motion without drift, and . This is a Gaussian process with correlogram • Account for nonstationarity by letting the kernel b vary with location:

40. Kernel averaging • Fuentes (2000): Introduce orthogonal local stationary processes Zk(s), k=1,...,K, defined on disjoint subregions Sk and construct • where wk(s) is a weight function related to dist(s,Sk). Then • A continuous version has

41. Simplifying assumptions in space-time models • Temporal stationarity • seasonality • decadal oscillations • Spatial stationarity • orographic effects • meteorological forcing • Separability • C(t,s)=C1(t)C2(s)

42. SARMAP revisited • Spatial correlation structure depends on hour of the day (non-separable):

43. Bruno’s seasonal nonseparability • Nonseparability generated by seasonally changing spatial term • Z1 large-scale feature • Z2 separable field of local features • (Bruno, 2004)

44. A non-separable class of stationary space-time covariance functions • Cressie & Huang (1999): • Fourier domain • Gneiting (2001): f is completely monotone if (-1)n f (n) ≥ 0for all n. Bernstein’s theorem : for some non-decreasing F. • Combine a completely monotone function and a function y with completely monotone derivative into a space-time covariance

45. A particular case a=1/2,g=1/2 a=1/2,g=1 a=1,g=1/2 a=1,g=1

46. Uses for surface estimation • Compliance • exposure assessment • measurement • Trend • Model assessment • comparing (deterministic) model to data • approximating model output • Health effects modeling

47. Health effects • Personal exposure (ambient and non-ambient) • Ambient exposure • outdoor time • infiltration • Outdoor concentration model for individual i at time t

48. Seattle health effects study • 2 years, 26 10-day sessions • A total of 167 subjects: • 56 COPD subjects • 40 CHD subjects • 38 healthy subjects • (over 65 years old, non-smokers) • 33 asthmatic kids • A total of 108 residences: • 55 private homes • 23 private apartments • 30 group homes

49. Ogawa sampler PUF HPEM pDR

50. T/RH logger CO2 monitor Ogawa sampler CAT Nephelometer HI Quiet Pump Box