Uniﬁed Geostatistical Modeling for Data Fusion and Spatial Heteroskedasticity with R Package ramps

This article illustrates usage of the ramps R package, which implements the reparameterized and marginalized posterior sampling (RAMPS) algorithm for complex Bayesian geostatistical models. The RAMPS methodology allows joint modeling of areal and point-source data arising from the same underlying spatial process. A reparametrization of variance parameters facilitates slice sampling based on simplexes, which can be useful in general when multiple variances are present. Prediction at arbitrary points can be made, which is critical in applications where maps are needed. Our implementation takes advantage of sparse matrix operations in the Matrix package and can provide substantial savings in computing time for large datasets. A user-friendly interface, similar to the nlme mixed eﬀects models package, enables users to analyze datasets with little programming eﬀort. Support is provided for numerous spatial and spatiotemporal correlation structures, user-deﬁned correlation structures, and non-spatial random eﬀects. The package features are illustrated via a synthetic dataset of spatially correlated observation distributed across the state of Iowa, USA.


Introduction
Spatial data, either areal or geostatistical (point-referenced), are becoming increasingly utilized in the study of many scientific fields due to the accessibility of data monitoring systems and associated datasets. When both types of data are available from the same underlying spatial process, computationally efficient and statistically sound methods are needed for their joint analysis. Markov chain Monte Carlo (MCMC) is a very powerful tool often used for the Bayesian analysis of such data. However, MCMC efficiency can be diminished by substantial autocorrelation in values of the model parameters sampled from the posterior distribution. Yan, Cowles, Wang, and Armstrong (2007) recently proposed a reparameterized and marginalized posterior sampling (RAMPS) algorithm which leads to lower autocorrelation in MCMC samples for Bayesian spatiotemporal geostatistical modeling. The RAMPS algorithm has been further extended to a unified framework of linear mixed models (Cowles, Yan, and Smith 2007) that allows fusion of data obtained at different resolutions (areal and point-referenced) and spatial heteroskedasticity. The general framework also covers cases where prediction at arbitrary sites and non-spatial random effects are needed. This article describes the implementation of the RAMPS algorithm in the R package ramps (Smith, Yan, and Cowles 2008) and illustrates its use with a synthetic dataset.
Existing R packages for geostatistical analysis include fields (Fields Development Team 2006), geoR (Ribeiro and Diggle 2001), geoRglm (Christensen and Ribeiro 2002), gstat (Pebesma and Wesseling 1998), sgeostat (Majure and Gebhardt 2007), spatial (Venables and Ripley 2002), and spBayes (Finley, Banerjee, and Carlin 2007). The fields, gstat, sgeostat, and spatial packages rely on frequentist kriging for modeling and prediction of geostatistical data. The geoR (and associated geoRglm package for generalized linear models) and spBayes packages offer routines to fit Bayesian geostatistical models. These packages do not accommodate combined analysis of point-source and areal data, which is one of the unique features of the ramps package. Furthermore, the spBayes package is not tailored to yield MCMC samples with lower auto-correlations, which may be critically important in analyzing large datasets. The geoR package attains independent posterior samples at the expense of discretizing the prior and posterior densities with respect to the spatial variance and correlation parameters.
The starting point for our unified geostatistical model is the basic RAMPS algorithm for point-source data only, described first in Yan et al. (2007). Consider geostatistical observations in a spatial domain D: {Y (s i ) : s i ∈ D, i = 1, . . . , n}, and let Y = {Y (s 1 ), . . . , Y (s n )} . A Gaussian geostatistical model for Y consists of spatial trend, spatial correlation, and measurement error: where β is a p × 1 vector of coefficients for covariate matrix X = {X (s 1 ), . . . , X (s n )} , Z is a n × 1 vector capturing the spatial correlation, and ε is a n × 1 vector of independent and identically distributed measurement errors. The distribution of Z is multivariate normal with mean zero and covariance matrix σ 2 z R(φ), where R(φ) is the correlation matrix as a function of parameter vector φ.
The RAMPS algorithm of Yan et al. (2007) includes two steps -reparameterization and marginalization -before drawing samples from the posterior density. The reparameterization step rewrites the model as where σ 2 = σ 2 z + σ 2 e and κ = σ 2 e /σ 2 . Letting θ = (φ, κ, σ 2 , β), the marginalization step factors the posterior density p(θ|Y ) as With appropriate prior distributions for elements in θ, the second and third distributions on the right hand side can be shown to be inverse gamma and Gaussian, respectively, which makes sampling from them very easy. Conversely, the first distribution is difficult to sample from, and Yan et al. (2007) used slice sampling for this critical step. Cowles et al. (2007) subsequently generalized the basic RAMPS algorithm to accommodate the following data complexities and research needs: 1. Fusion of areal data and point-referenced data in a single model. Such data fusion combines data from different sources and resolutions to make more precise statistical inferences, which oftentimes is in terms of narrower credible sets for parameter estimation and prediction.
2. Non-spatial random effects. An example of such non-spatial random effects is the radon data analysis of Smith and Cowles (2007), where many sites have multiple measurements and a site-specific random effect is needed.
3. Multiple variances for each variation source. In fact, data fusion naturally demands multiple variances in the measurement error process for different data sources. The general model framework not only meets this demands but also allows the underlying spatial process to have different variances at different locations; that is, spatial heteroskedasticity.
4. Prediction at arbitrary sites, measured or unmeasured. The RAMPS algorithm can be carried out with very little change in formulation using the method of structured hierarchical models (Hodges 1998;Sargent, Hodges, and Carlin 2000).
All these capabilities are implemented in the ramps package, which is publicly available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=ramps.
The ramps package offers a complete set of tools for the conduct of Bayesian geostatistical analysis of large, complex spatial datasets using the RAMPS algorithm. Its unique features are summarized in Table 1 from the aspects of modeling, computing, correlation structures, and user-interface. Note that the correlation structures summarized in Table 1 are a superset of those provided in nlme (Pinheiro and Bates 2000) with additional support in ramps for great-circle distance as a distance metric option.
This article is organized as follows. Section 2 presents a unified geostatistical model framework that incorporates the aforementioned generalizations. Section 3 discusses some implementational details of the ramps package as well as its user interface. Section 4 illustrates the use of the package through a working example. Section 5 reports a performance evaluation of the package in comparison with packages spBayes and geoR in the context of fitting a simple geostatistical model. A discussion concludes in Section 6.

Unified geostatistical model
Let Y = (Y a , Y p ) be a concatenated vector of areal observations Y a = {Y a,1 , . . . , Y a,na } and point-referenced observations Y p = {Y p,1 , . . . , Y p,np } , where n a + n p = n is the total number of observations. The unified Gaussian geostatistical model is Modeling • Joint modeling of data from multiple sources (point-source, areal, or both).
• Non-spatial random effects as in a linear mixed model. • Multiple variances for each variation source (measurement error, spatial, and random effects). • Prediction at measured or unmeasured sites.

Computing
• Efficient MCMC sampling with the RAMPS algorithm.
• Sparse matrix operation exploited for large datasets.
User interface • Easy-to-use model specification.
• Three-dimensional spatial plotting of results. where X, W , and K are design matrices for fixed effects β (p × 1), non-spatial random effects γ (q × 1), and spatial random effects Z (S × 1), respectively. The matrix K is defined by otherwise.
In the case of a point-referenced observation, one measurement site contributes to Y i , and thus N i = 1. Conversely, multiple measurement sites contribute to an areal observation Y i . The model defines such an observation as the average over N i > 1 sites. If the actual numbers or locations of contributing sites are unknown, then a uniform grid of spatial sites may be used as an approximation. Accordingly, the N i will be roughly proportional to the area of the region over which Y i is an average. The finer the grid of sites; the closer the proportionality will be. In summary, the model formulation (2) accommodates point-referenced data, areal data, multiple measurements, and non-spatial random effects.
Data fusion is made possible in model (2) through the allowance of both areal and pointreferenced data. When both types are included simultaneously, a common underlying spatial process Z(s) is assumed, and the design matrix K maps Z contributions to the observed data. For point-referenced data at site s, the contribution from Z is simply Z(s). For data averaged over an area A, the contributions are from {Z(s) : s ∈ G, s ∈ A}, where G is a grid of sites defined over the region from which the data are collected. In practice, the spatial random effects Z in model (2) contain realizations of the spatial process Z(s) at all unique pointreferenced and grid sites. The fineness of the grid G can be tuned, depending on the scientific question and computational resources available. Note that Z can also contain realizations at sites that do not contribute to any observed data but at which prediction is of scientific interest, in which case, the corresponding rows in K will consist of zeros; a formulation for this purpose will be presented at the end of this section.
Heteroskedasticity is accommodated by allowing variances to differ across the non-spatial random effects, spatial measurement sites, and individual measurement types. Suppose that there are L γ different variances for the non-spatial random effects σ 2 γ,i , i = 1, . . . , L γ ; L Z spatial variances σ 2 Z,i , i = 1, . . . , L Z ; and L measurement error variances σ 2 ,i , i = 1, . . . , L . Further, let r k , k = 1, . . . , q, be an integer between 1 and L γ indicating the corresponding random effects variance for γ k . Likewise, let v j , j = 1, . . . , S, indicate the spatial variance for observations from site j, and m i , i = 1, . . . , n the measurement error variance for observation Y i . We construct vectors for componentwise variances of γ, Z, and , respectively, as . . , n, is a weight associated with observation i. In the ramps package, users may specify the weighting values or accept the default values of 1 for point-source and N i for areal observations. Assuming exchangeability of random effects, we have Σ is a spatial correlation matrix with parameter vector φ. This setup is general and allows modeling for spatiotemporal data.
The variance parameters are reparameterized to facilitate the marginalization of the posterior density in the RAMPS algorithm. Concatenate the vectors of measurement error variances, spatial variances, and random effects variances for a total of F = L γ + L Z + L variance parameters, σ 2 1 , . . . , σ 2 F . If there is one measurement-error variance, one spatial variance, and no random effects variances, then σ 2 1 ≡ σ 2 e and σ 2 2 ≡ σ 2 z as in the special case of Yan et al. (2007). Our reparameterization is in terms of κ = {κ 1 , . . . , κ F } and σ 2 tot , where Note that κ F ≡ 1 − F −1 j=1 κ j and is not a free parameter to be estimated. Let and Ω = diag(κ ). Then the likelihood can be specified as where Ω = W Ω γ W + KΩ Z K + Ω . Cowles et al. (2007) derived the factors of the posterior density p(φ, κ|Y ), p(σ 2 tot |φ, κ, Y ), and p(β|φ, κ, σ 2 tot , Y ). The prior distributions are inverse gamma on σ 2 j , j = 1, . . . , F , flat on β, and uniform for φ with appropriately chosen bounds. A challenge presented in sampling from p(φ, κ|Y ) is the constraint that κ has support on the standard (F − 1)-simplex The RAMPS procedure can be modified to accommodate prediction at arbitrary sites. Partition Z as (Z p , Z u ) , where Z p is a vector spatial random effects at sites for which prediction is desired, and Z u is a vector of spatial random effects at sites for which no prediction is desired. Sampling of β and Z p can be done in a batch by partitioning and rearranging the matrix K such that KZ = (K p , K u )(Z p , Z u ) . Similar to Sargent et al. (2000), a stacked linear model can be obtained as This extension simply revises the likelihood expression in equation (3) as and the RAMPS algorithm can be applied to sample B = (β , Z p ) .

Correlation structures
Characteristic of geostatistical models is the specification of correlation as a function of the distance between, and possibly orientation of, geographic points in the spatial domain. Our model as implemented in the ramps package allows spatial correlation of the general form where c (s i , s i ; φ) is a function of the distance between sites s i and s i and the parameter vector φ. We provide metrics for the calculation of spatial distance as great-circle distance, Euclidean distance, maximum distance, and sum of absolute differences. Available parametric correlation functions are summarized in Table 3.1. Usage is consistent across the correlation functions, and spatial covariates, such as longitude and latitude, are allowed in the formula specification; see Section 4 for illustration. These are extensions of the nlme spatial correlation structures and offer users a consistent interface for geostatistical model specification in the ramps package. The spatial correlation structures in nlme are not directly used because they do not allow great circle distance, which is very commonly needed for spatial data.
In addition to the supplied functions, users can create their own correlation structures for use with the package by defining a new corSpatial class and accompanying constructor, corMatrix, and coef method functions. Examples can be found in the ramps source code.

Model fitting interface
The main user-level function for geostatistical model fitting in the ramps package is georamps whose definition is given below.

Spatiotemporal correlation corRExp2
Exponential corRExpwr2 Powered exponential corRExpwr2Dt Temporally-integrated powered exponential This function implements the RAMPS algorithm for generating samples from the posterior distribution of the model parameters in geostatistical model (2). Model specification of the fixed effects (fixed), random effects (random), and spatial correlation (correlation) arguments mirrors those in package nlme. Data fusion and heteroskedasticity are specified by two separate arguments described as follows.
The optional argument aggregate allows specification of grid sites over which areal observations are assumed to be aggregated. It is fed by a list of elements: grid -a data frame of coordinates to use for Monte Carlo integration over geographic blocks at which areal measurements are available; and blockid -a character string specifying the column by which to merge the areal measurements in the data (data) with the grid coordinates in grid. Merging is performed only for blockid values that are common to both datasets. All observations in data are treated as point-source measurements by default.
The argument variance specifies the types of the multiple variances for each variation source. It is supplied an optional list of one-sided formulas, each of the form~g where g defines the grouping factors for: fixed -measurement error variances V ; random -random effects error variances V γ ; and spatial -spatial variances V Z . A single variance is assumed in each case by default.
Another important argument is control, which controls the fitting process through initial parameter values and prior distributions. Its value must be a ramps.control object generated from the ramps.control function described next.

Fitting control
The ramps.control function is defined as thus: ramps.control(iter = 1000, beta, sigma2.e, phi, sigma2.z, sigma2.re, z.monitor = TRUE, mpdfun = c("mpdbeta", "mpdbetaz"), file) It collects from the user the number of desired MCMC iterations (iter), the prior distribution for model parameters (see below), sites at which prediction is needed (z.monitor), and file names (file) for the output of monitored sample values to external files.
Prior distributions need to be specified with the param function for all model parametersfixed effects beta, spatial correlation parameter phi, variance parameters for measurement errors sigma2.e, spatial random effects sigma2.z, and non-spatial random effects sigma2.re.
Hyperparameters of the prior distributions are passed through the '...' argument.
Tuning of the initial hypercube/simplex sizes for slice sampling is controlled via the tuning argument. For spatial correlation parameters φ, the tuning parameters can take on positive values and define the factors by which the associated uniform priors are multiplied to define the initial hyperrectangle size. Tuning parameter values to control the initial simplex size for κ may be given in the variance specifications and must be between 0 and 1. The smallest tuning parameter among the variance specifications will be used in tuning the slice algorithm for κ. Smaller values of the tuning parameters will lead to faster sampling at the expense of higher autocorrelation in sampler output.

Model comparisons
Bayesian model fits are commonly compared with the deviance information criterion (DIC), defined as where ϑ denotes the model parameters, D (Y | ϑ) the associated deviance function, and p D a measure of the effective number of model parameters (Spiegelhalter, Best, Carlin, and van der Linde 2002). Smaller values of the deviance function D are indicative of models that provide better fits to the data. Since this function necessarily decreases as the number of parameters increases, the effective number of parameters p D is included in the criterion as a penalty term. Consequently, comparisons based on the DIC aim to strike a balance between model goodness-of-fit and parsimony. Preference is given to models with smaller DIC values.
In the context of our geostatistical model, deviance is calculated as a function of the data likelihood Y ∼ N (XB, σ 2 tot Ω).
Hence, DIC results are dependent upon the form of the likelihood function used. Two different likelihood formulations are implemented for DIC calculation in the ramps package. In the first, X = X, B = β, and Ω = W Ω γ W + KΩ Z K + Ω so that the random effects are integrated out of the mean structure, leaving only the fixed covariates and associated β parameters. In the second formulation, both fixed and spatial random effects are included in an SMCMC specification of the likelihood, such that The mpdfun argument of the ramps.control function determines which one of the two likelihood formulations is used for DIC calculations -values of "mpdbeta" and "mpdbetaz" correspond to the first and second formulations, respectively. Some guidance on choosing between the two options is provided by the simulation results of Section 4.

Computational performance
MCMC algorithms for geostatistical models are computationally intensive because the spatial correlation matrix must be updated and decomposed at each sampler iteration. The latter, more expensive operation is of computational order equal to the number of unique spatial sites cubed. Thus, the speed of MCMC algorithms is highly dependent on the performance of the linear algebra routines used. The ramps package relies upon core R functions for linear algebra operations. These core function interface with routines provided by the BLAS and LAPACK linear algebra libraries against which R is linked at runtime. Consequently, one effective way to improve substantially runtime performance of the MCMC algorithm in ramps is to link R against optimized versions of BLAS and LAPACK, such as those provided by ATLAS (Whaley, Petitet, and Dongarra 2001), ACML (Advanced Micro Devices, Inc. 2008), Goto BLAS (Goto and van de Geijn 2006), or Intel MKL (Intel 2008). Details on the configuration of R for these alternative linkages are provided in R's accompanying "Installation and Administration" guide.
Other optimization techniques have been applied to the code directly to improve performance. For example, the K, W , X and Ω matrices in the likelihood tend be very sparse, hence suggesting the use of sparse matrix libraries as one way to accelerate computations. Recent versions of the Matrix package (Bates and Maechler 2007) provide interfaces to the sparse matrix libraries of Davis (2006) and are used in the implementation of our ramps package. As sample size increases, the speed advantages of the Matrix routines for sparse matrix operations are well worth the implementation effort. Another example is the choice of marginalized posterior density p(φ, κ|Y ) to be used in the SIMPLICE step of the MCMC algorithm. With the mpdfun argument of the ramps.control function, users can specify that the density be marginalized with respect to the β parameters only ("mpdbeta") or with respect to both the β and Z parameters ("mpdbetaz"). The latter choice is recommended as a way to increase algorithm speed when there are multiple observations per measurement site. Regardless of the marginalization choice, samples are drawn from the same posterior distribution.

Working example
To illustrate the use of ramps for joint analysis of areal and point-source observations, a synthetic dataset was generated from model (2) using the county structure in the state of Iowa, USA. There are n a = 99 counties in the state. Areal observations were generated as county averages from a uniform grid of 391 sites -approximately 4 sites per county. In addition, we generated n p = 600 point-source observations from a set of n s = 300 unique sites sampled from a uniform distribution in Iowa. The multiple point-source measurements per site allow for site-specific non-spatial random effects in the analysis. In Figure 1 the grid of 391 sites is depicted with circles and the 300 point-source measurement sites with dots.
An underlying spatial process was generated from a multivariate normal distribution using an exponential correlation structure with range parameter φ = 10 and variance σ 2 z = 0.36. A value of φ = 10 implies that the correlation between two sites drops to 0.05 at a distance of about 30 miles. Measurement errors were generated with variances σ 2 ε,0 = 0.25 for point-source data and σ 2 ε,1 = 0.09 for areal data. Site-specific non-spatial random effects were generated with variance σ 2 γ = 0.16. One fixed effects covariate areal is included as an indicator for areal observations. Its β coefficient was set equal to 0.5.
Simulated data are stored in the data frame simIowa, with columns y for the areal and point-source observations, areal, lon and lat giving the longitude and latitude coordinates, siteId as a unique site identifier, and weights containing weighting values for measurement error variances. In order to combine the two types of observations in one dataset, lon, lat, and siteId are assigned missing NA values for areal observations. A separate grid of measurements sites for areal observations must be supplied to the georamps function. The latitude and longitude coordinates of the 391 uniform grid sites in our example are stored in the data frame simGrid as variables lon and lat. An additional indexing variable id is included in both simGrid and simIowa for the purpose of matching grid sites to their respective areal observations. The code below creates a control object of model fitting parameters that must be supplied to the georamps function. Prior distributions on θ are: Unif(1, 60) for φ, IG(0.01, 0.01) for σ 2 ,1 , σ 2 ,2 , σ 2 z , and σ 2 γ , and flat for β. Also specified are the number of MCMC iterations (iter), coordinates of sites at which prediction is desired (z.monitor), and optional names of external files to which to save MCMC output for model parameters and spatial random effects (file).
The joint analysis of areal and point referenced data can now be performed with a call to georamps: R> fit.fusion <-georamps(fixed = y~areal, + random =~1 | siteId, + correlation = corRExp(form =~lon + lat, metric = "haversine"), + variance = list(fixed =~areal), weights = weights, + data = simIowa, + aggregate = list(grid = simGrid, blockid = "id"), + control = control.fusion) The model has one covariate (areal) as a fixed effect and a site-specific (siteId) random effect. The spatial correlation structure is exponential, corRExp, with spatial distance computed as great-circle distance (haversine). Of note is that, when the haversine metric is used, the order of variables must be such that longitude is first and latitude is second. The argument variance specifies grouping factors for each variance component associated with the measurement errors (fixed), non-spatial random effects (random), and spatial random effects (spatial). The argument aggregate is simply a list which gives the grid from which the areal data are assumed to be obtained and the name of the variable with which the grid and observed data can be merged. The aggregate argument would be omitted if analyzing point-source-only data.
We ran 1,100 MCMC iterations and discarded the first 100 as a burn-in sequence. The remaining 1,000 iterations were used for posterior inference. For instance, posterior summaries for the joint analysis were obtained with the code given below.
R> fit.fusion1000 <-window(fit.fusion, iter = 101:1100) R> summary(fit.fusion1000) For comparison, we performed analyses separately for both the point-source only data and the areal only data. The code can be written by modifying that given previously for the fused data analysis and is illustrated in the package help files. Table 4 summarizes the percentiles of the posterior samples from the three analyses. Results in Table 4 indicate that the joint analysis gives narrower 95% credible intervals for parameters common to all analyses; e.g., φ and σ 2 Z . Additionally, we used the ramps DIC function to evaluate model fits under different specifications of the spatial correlation structure.

R> DIC(fit.fusion1000)
Results are presented in Table 4 for the two implemented forms of the data likelihoodmpdbeta and mpdbetaz. Recall that smaller DIC values are suggestive of more desirable models. For selecting among correlation structures, better discrimination is provided by the mpdbetaz likelihood. In particular, preference is shown here for the exponential structure which was used to simulate the data. Finally, we excluded the fixed effects indicator for areal observations from the joint model, obtaining DIC values of 1488.5 (mpdbeta) and 1921.3 (mpdbetaz); both confirming the importance of the indicator in the analysis.   Table 4: Deviance information criterion values for different spatial correlation specifications, with effect number of parameters given in parentheses. Values are based on likelihood formulations with only fixed effects (mpdbeta) and both fixed and spatial random effects (mpdbetaz) in the mean structure. The numerical approximation of Gneiting (1999) is used for the Gaussian correlation.
Mapping of the spatial distribution is often of particular interest. There are two ways to get MCMC samples of spatial random effects. The first way is set z.monitor in function ramps.control equal to TRUE or a data frame of coordinates at which prediction is needed. This way is designed for sites that contribute to the observed data. The second way is to use the predict method on the ramps object returned by georamps. This way is designed for sites that do not contributed to the observed data and is particularly useful when prediction on a new grid of sites is needed to draw maps after analyses of point-source data. For example, consider the new grid of sites created with the code given below.
In addition to color image maps of the spatial distribution, the plot function provides a type argument that allows for the construction of wireframe ("w") and contour ("c") maps, as shown in the code below and corresponding Figure 3. Note that the plot function easily accommodates other choices of color palettes, such as the one derived here from the HCL color space and described by Zeileis, Hornik, and Murrell (2007).
MCMC samplers were run for 1,000 iterations using each package, starting from the same initial values of all parameters. All three packages give very similar quantiles of the MCMC samples. The autocorrelations, however, are quite different, which is reflected in terms of "effective sample size" (ESS) (Kass, Carlin, Gelman, and Neal 1998), an established metric for comparing performance of MCMC algorithms. ESS is the number of independent samples that would carry the same amount of information as the available correlated samples. For a given number of MCMC sampler iterations, the higher the autocorrelation in sampler output for a particular parameter, the smaller the resulting effective sample size. Speed of sampling algorithms can be compared fairly in terms of the effective samples per unit run time. The effectiveSize function in the R package coda (Plummer, Best, Cowles, and Vines 2006) was used to calculate the effective sample size for each parameter from the output of 1,000 MCMC iterations generated with each package.
The ESS and ESS per minute for the 1,000 MCMC samples using the three packages are summarized in Table 5. For the regression coefficient β, all three packages do well, giving 1,000 (or 929 for spBayes) independent draws. For the spatial variance σ 2 z , the most difficult parameter to estimate, the ramps packages produces a sample worth 272.8 independent draws, about 8 times as many as the spBayes package gives (36.1). When time is taken into consideration, the ramps package takes 73.4 minutes while the spBayes packages takes 52.8 minutes. The ramps package gets 5.5 times as many ESS per minute as the spBayes package does. The geoR package produces higher ESS and ESS/min than the ramps package, but recall that the posterior samples are obtained on a grid. The distribution is discrete and the posterior density is jagged.

Discussion
The ramps package enables Bayesian geostatistical analysis with the very general class of models described by (2). As exemplified in the performance evaluation, its implementation of the RAMPS algorithm provides the advantage of low autocorrelation in MCMC output and therefore more effective samples per unit time than competing methods. The SIMPLICE algorithm which performs slice sampling based on simplexes can be generally useful for cases where multiple variances are present (He, Hodges, and Carlin 2007). As a geostatistical tool, the package also provides smooth maps for either point-source or areal observations. Furthermore, users have full control over specification of the grid from which areal observation are assumed to be drawn.
In our experiments, the spatial correlation parameter φ has usually been hardest to estimate and its posterior sample autocorrelation highest among all parameters. Conversely, the fixed effects are easiest to estimate and their posterior samples almost independent. This observation shows the importance of tuning the size of the initial hyperrectangle for φ and simplex for κ in the slice sampling by setting tuning in param when defining the control object with ramps.control. For large datasets, one may wish to choose a larger tuning parameter for φ and a smaller tuning parameter for κ such that sampling of φ traverses the sample space more quickly.
The efficiency of our RAMPS algorithm is determined by the autocorrelation in sampling the marginalized density p(φ, κ|Y ). For lower dimensional (φ, κ), it is possible to evaluate the density on a grid first, which can then be used to produce close-to-independent samples, as is done in geoR. However, for higher dimensional (φ, κ), such as arise with spatiotemporal correlation functions and multiple sources of variance, discrete grid evaluation may not be feasible. The RAMPS algorithm and package are designed to accommodate the latter, more general setting. Nevertheless, an adaptive MCMC procedure which takes advantage of past marginalized density evaluations is worth future investigation.