Douglas Nychka

Semi-nonparametric Maximum Likelihood Estimation

Often maximum likelihood is the method of choice for fitting an econometric model to data but cannot be used because the correct specific ation of (multivariate) density that defines the likelihood is unknown. In this situation, simply put the density equal to a Hermite series and apply standard finite dimensional maximum likelihood methods. Model parameters and nearly all aspects of the unknown density itself will be estimated consistently provided that the length of the series increases with sample size. The rule for increasing series length can be data dependent. The method is applied to nonlinear regression with sample selection. Copyright 1987 by The Econometric Society.

Covariance Tapering for Interpolation of Large Spatial Datasets

Interpolation of a spatially correlated random process is used in many scientific areas. The best unbiased linear predictor, often called a kriging predictor in geostatistical science, requires the solution of a (possibly large) linear system based on the covariance matrix of the observations. In this article, we show that tapering the correct covariance matrix with an appropriate compactly supported positive definite function reduces the computational burden significantly and still leads to an asymptotically optimal mean squared error. The effect of tapering is to create a sparse approximate linear system that can then be solved using sparse matrix algorithms. Monte Carlo simulations support the theoretical results. An application to a large climatological precipitation dataset is presented as a concrete and practical illustration.

Bayesian Spatial Modeling of Extreme Precipitation Return Levels

Quantification of precipitation extremes is important for flood planning purposes, and a common measure of extreme events is the r-year return level. We present a method for producing maps of precipitation return levels and uncertainty measures and apply it to a region in Colorado. Separate hierarchical models are constructed for the intensity and the frequency of extreme precipitation events. For intensity, we model daily precipitation above a high threshold at 56 weather stations with the generalized Pareto distribution. For frequency, we model the number of exceedances at the stations as binomial random variables. Both models assume that the regional extreme precipitation is driven by a latent spatial process characterized by geographical and climatological covariates. Effects not fully described by the covariates are captured by spatial structure in the hierarchies. Spatial methods were improved by working in a space with climatological coordinates. Inference is provided by a Markov chain Monte Carlo algorithm and spatial interpolation method, which provide a natural method for estimating uncertainty.

Bayesian Confidence Intervals for Smoothing Splines

Abstract The frequency properties of Wahbas Bayesian confidence intervals for smoothing splines are investigated by a large-sample approximation and by a simulation study. When the coverage probabilities for these pointwise confidence intervals are averaged across the observation points, the average coverage probability (ACP) should be close to the nominal level. From a frequency point of view, this agreement occurs because the average posterior variance for the spline is similar to a consistent estimate of the average squared error and because the average squared bias is a modest fraction of the total average squared error. These properties are independent of the Bayesian assumptions used to derive this confidence procedure, and they explain why the ACP is accurate for functions that are much smoother than the sample paths prescribed by the prior. This analysis accounts for the choice of the smoothing parameter (bandwidth) using cross-validation. In the case of natural splines an adaptive method for avo...

Climate spectra and detecting climate change

Part of the debate over possible climate changes centers on the possibility that the changes observed over the previous century are natural in origin. This raises the question of how large a change could be expected as a result of natural variability. If the climate measurement of interest is modelled as a stationary (or related) Gaussian time series, this question can be answered in terms of (a) the way in which change is estimated, and (b) the spectrum of the time series. These computations are illustrated for 128 years of global temperature data using some simple measures of change and for a variety of possible temperature spectra. The results highlight the time scales on which it is important to know the magnitude of natural variability. The uncertainties in estimates of trend are most sensitive to fluctuations in the temperature series with periods from approximately 50 to 500 years. For some of the temperature spectra, it was found that the standard error of the least squares trend estimate was 3 times the standard error derived under the naïve assumption that the temperature series was uncorrelated. The observed trend differs from zero by more than 3 times the largest of the calculated standard errors, however, and is therefore highly significant.

Design of Air-Quality Monitoring Networks

Where should ozone be measured? It is well accepted that high levels of ozone are not only damaging to human health but also reduce crop yield and damage vegetation.1 However, the continuous measurement of ozone at a location is relatively expensive and so the number and locations of instruments need to be chosen judiciously. This question, although deceptively simple, raises a host of fundamental issues. Most importantly, how can we infer ozone levels at places where measurements are not made? What does it mean to measure ozone well and how many monitoring instruments are really necessary?

Estimating the Lyapunov Exponent of a Chaotic System with Nonparametric Regression

Abstract We discuss procedures based on nonparametric regression for estimating the dominant Lyapunov Exponent λ1 from time series data generated by a nonlinear autoregressive system with additive noise. For systems with bounded fluctuations, λ1 > 0 is the defining feature of chaos. Thus our procedures can be used to examine time series data for evidence of chaotic dynamics. We show that a consistent estimator of the partial derivatives of the autoregression function can be used to obtain a consistent estimator of λ1. The rate of convergence we establish is quite slow; a better rate of convergence is derived heuristically and supported by simulations. Simulation results from several implementations—one “local” (thin-plate splines) and three “global” (neural nets, radial basis functions, and projection pursuit)—are presented for two deterministic chaotic systems. Local splines and neural nets yield accurate estimates of the Lyapunov exponent; however, the spline method is sensitive to the choice of the emb...

Noise and Nonlinearity in Measles Epidemics: Combining Mechanistic and Statistical Approaches to Population Modeling

We present and evaluate an approach to analyzing population dynamics data using semimechanistic models. These models incorporate reliable information on population structure and underlying dynamic mechanisms but use nonparametric surface‐fitting methods to avoid unsupported assumptions about the precise form of rate equations. Using historical data on measles epidemics as a case study, we show how this approach can lead to better forecasts, better characterizations of the dynamics, and a better understanding of the factors causing complex population dynamics relative to either mechanistic models or purely descriptive statistical time‐series models. The semimechanistic models are found to have better forecasting accuracy than either of the model types used in previous analyses when tested on data not used to fit the models. The dynamics are characterized as being both nonlinear and noisy, and the global dynamics are clustered very tightly near the border of stability (dominant Lyapunov exponent λ ≈ 0). However, locally in state space the dynamics oscillate between strong short‐term stability and strong short‐term chaos (i.e., between negative and positive local Lyapunov exponents). There is statistically significant evidence for short‐term chaos in all data sets examined. Thus the nonlinearity in these systems is characterized by the variance over state space in local measures of chaos versus stability rather than a single summary measure of the overall dynamics as either chaotic or nonchaotic.

Multivariate Bayesian analysis of atmosphere-ocean general circulation models

Numerical experiments based on atmosphere–ocean general circulation models (AOGCMs) are one of the primary tools in deriving projections for future climate change. Although each AOGCM has the same underlying partial differential equations modeling large scale effects, they have different small scale parameterizations and different discretizations to solve the equations, resulting in different climate projections. This motivates climate projections synthesized from results of several AOGCMs’ output. We combine present day observations, present day and future climate projections in a single highdimensional hierarchical Bayes model. The challenging aspect is the modeling of the spatial processes on the sphere, the number of parameters and the amount of data involved. We pursue a Bayesian hierarchical model that separates the spatial response into a large scale climate change signal and an isotropic process representing small scale variability among AOGCMs. Samples from the posterior distributions are obtained with computer-intensive MCMC simulations. The novelty of our approach is that we use gridded, high resolution data covering the entire sphere within a spatial hierarchical framework. The primary data source is provided by the Coupled Model Intercomparison Project (CMIP) and consists of 9 AOGCMs on a 2.8 by 2.8 degree grid under several different emission scenarios. In this article we consider mean seasonal surface temperature and precipitation as climate variables. Extensions to our model are also discussed.

A multi-resolution Gaussian process model for the analysis of large spatial data sets

We develop a multiresolution model to predict two-dimensional spatial fields based on irregularly spaced observations. The radial basis functions at each level of resolution are constructed using a Wendland compactly supported correlation function with the nodes arranged on a rectangular grid. The grid at each finer level increases by a factor of two and the basis functions are scaled to have a constant overlap. The coefficients associated with the basis functions at each level of resolution are distributed according to a Gaussian Markov random field (GMRF) and take advantage of the fact that the basis is organized as a lattice. Several numerical examples and analytical results establish that this scheme gives a good approximation to standard covariance functions such as the Matérn and also has flexibility to fit more complicated shapes. The other important feature of this model is that it can be applied to statistical inference for large spatial datasets because key matrices in the computations are sparse. The computational efficiency applies to both the evaluation of the likelihood and spatial predictions.