The motivation for this dissertation is the need that often arises in spatial settings to perform a data transformation to achieve a stationary process and/or variance stabilization. The transformation may be a non-linear transformation, and the desired predictand may be multivariate in that it is necessary to interpolate predictions at multiple sites. We assume the underlying spatial model is a Gaussian random field with a parametrically specified covariance structure, but that the predictions of interest are for multivariate nonlinear functions of the Gaussian field. This induces new complications in the spatial interpolation known as kriging. For instance, it is no longer possible to derive the predictive distribution function in closed form. The underlying process of a spatial model is a stochastic process, and spatial prediction techniques are founded on the assumption that the realizations come from a Gaussian process with a known covariance structure, and known, specified parameters. A difficulty that arises with traditional kriging methods is the fact that the standard formula for the mean squared prediction error does not take into account the estimation of the covariance parameters. This generally leads to underestimated prediction errors, even if the model is correct. Smith and Zhu (2004) establish a second-order expansion for predictive distributions in Gaussian processes with estimated covariances. Here, we establish a similar expansion for multivariate kriging with non-linear predictands. Bayesian methods provide a possible resolution to errors encountered through employing frequentist estimation techniques for obtaining spatial parameters. Bayesian analysis seeks to utilize prior and nonexperimental sources of information. Bayesian methods evaluate procedures for a repeated sampling experiment of parameters drawn from the posterior distribution given a set of observed data. An important property of Bayesian methods is the ability to deal with the uncertainty in a particular model. A Bayesian paradigm enables a more realistic assessment of the variability inherent in estimating parameters of interest. Here we explore a Laplace approximation to Bayesian techniques that provides an alternative to common iterative Bayesian methods, such as Markov Chain Monte Carlo. The theoretical Bayesian coverage probability bias for prediction intervals is computed and compared with the plug-in method using the restricted maximum likelihood estimates of the covariance parameters. The form of the Bayesian predictive distribution can be expressed as partial derivatives of the restricted log-likelihood. This leads to possible analytical evaluation, and a bootstrap method is explored to obtain predictions for general, non-linear predictands. The main results are asymptotic formulae for a general, non-linear predictand for the expected length of a Bayesian prediction interval, which has possible applications in network design, and for the coverage probability bias, which can lead to the development of a matching prior. The matching prior may be difficult to compute in practice. As an alternative, an asymptotic estimator is developed.