Simulating data for illustration
To illustrate our method, we simulate a small dataset. First, consider a unit square and randomly select observation locations. (Set spatial.dim=1 to consider data on the one-dimensional unit interval.)
Next, we define the covariance function of the field as well as the scale of the measurement error
We are now ready to simulate the field and visualize it as a sanity check.
We also create a grid of \(n.p\) locations at which we would like to make predictions.
More details on likelihood evaluation
Let’s take a closer look at how the likelihood is evaluated using Vecchia. Most importantly, we can specify a parameter, \(m\). Its value determines the number of “neighbours” of each point, or, in other words, how many other points a given point conditions on. The larger this parameter, the more accurate and expensive the approximation will be.
Note that the function vecchia_specify determines the general properties of the approximation, but it does not depend on the data or the specific parameter values. Hence, it does not have to be re-run when searching over different parameter values in an estimation procedure.
We can also compare the results to the exact likelihood:
In this case the approximation is very good. In general, \(m=20\) is a good value, and \(m\) should usually be between 10 and 40. For one-dimensional space, we can get good approximations even with \(m=5\) or smaller.
More details on spatial prediction
Similar to the previous section we next specify the approximation and indicate at which locations prediction is desired.
It is also possible to print the entire predictive covariance matrix. We do it here only for the purpose of illustration. If \(n.p\) is very large, this matrix might use up a lot of memory and we generally do not recommend plotting it directly.
Linear combinations
We might sometimes be interested in a linear combination of the predicted values. In particular, we can limit our attention to only a subset of our predictions. This can be accomplished by specifying the linear combination coefficients as a matrix. As an example, we assume we are only interested in predictions at the unobserved prediction locations (not at the first n observed locations):
As another example, we consider the overall mean of the process at all prediction locations. Using the vecchia_lincomb() function enables us to get the variance estimates easily.
Other GP approximations as special cases
By specifying appropriate options in vecchia_specify, we can do everything described above for several other GP approximations: Modified predictive process, FSA, MRA, latent, standard Vecchia
Setting \(M=1\) results in block full-scale approximation, specifically one with \(r_0 = \frac{m}{2}\) knots spread over the entire domain and the remaining locations being partitioned into blocks of size \(<\frac{m}{2}+1\).
Another popular existing approximation method, modified predictive process (MPP), can also be obtained by specifying appropriate parameter settings:
As we can see, MPP can be viewed as a special case of the multi-resolution approximation (MRA).
A general MRA is obtained my specifying all of its three parameters
mra.options.mra = list(r=c(10, 5, 5), M=2, J=2)
MRA_rJM = vecchia_specify(locs, m, 'maxmin', conditioning='mra', mra.options=mra.options.mra, verbose=TRUE)
#> Warning in get.mra.params(n, mra.options, m): M, r set for MRA. If parameter m
#> was given, it will be overridden
#> MRA params: m=22; J=2,2; r=10,5,7; M=2
We should note two things to note about this full specifiction of an MRA. First, providing all three \(r\),\(J\) and \(M\) overrides whatever value of \(m\) was provided. Second, in order to be able to place a knot at each point of the grid, the provided parameters might need to be adjusted.
Finally, we can also use the GPvecchia package to specify a Nearest Neighbour Gaussian Process (NNGP) approximation. This can be accomplished as shown below.
We can now easily compare different approximation methods and compare it with SGV and exact likelihood.
vecchia_likelihood(z,blockFS,covparms,nuggets)
#> [1] -51.62637
vecchia_likelihood(z,MPP,covparms,nuggets)
#> [1] -55.6042
vecchia_likelihood(z,MRA_rJM,covparms,nuggets)
#> [1] -51.30935
vecchia_likelihood(z,NNGP,covparms,nuggets)
#> [1] -51.41268
vecchia_likelihood(z, vecchia_specify(locs, m), covparms, nuggets)
#> [1] -51.40362
dmvnorm(z,mean=rep(0,n),sigma=Om0,log=TRUE)
#> [1] -51.441
Non-Gaussian data
Here we demonstrate how GPVecchia can fit a latent model to non-Gaussian data using the Vecchia-Laplace method. We simulate data by first generating a correlated latent field without noise, assuming the same covariance and locations generated earlier:
Then we sample a single non-Gaussian value for each latent value. The variability introduced by the sampling induces heteroskedasticity, in contrast the the constant noise added to the Gaussian case. Below we use a logistic model for binary data, but there are implementations for count and continuous positive (right-skewed) data as well.
data.model = "logistic"
# simulate data
if(data.model=='poisson'){
z = rpois(n, exp(y))
} else if(data.model=='logistic'){
z = rbinom(n,1,prob = exp(y)/(1+exp(y)))
} else if(data.model=='gamma'){
z = rgamma(n, shape = default_lh_params$alpha, rate = default_lh_params$alpha*exp(-y))
}else{
print('Error: Distribution not implemented yet.')
}
# plot simulated data, 1 or 2D
defpar = par(mfrow=c(1,2))
if(spatial.dim==1) {
plot(locs,y, main = "latent")
plot(locs,z, main = "observed")
} else {
fields::quilt.plot(locs,y, main = "Latent")
fields::quilt.plot(locs,z, main = "Observed")
}
Given the simulated data, we now can efficiently estimate the latent field by specifying the number of conditioning points \(m\) described earlier. Interweaved ordering is best for 1D data while response-first (‘zy’) ordering is best for higher dimensions.
With the approximated covariance structure, we can calculate the posterior estimate for the latent field using the Vecchia-Laplace method and plot the result. Pure Laplace approximation is included for comparison; even with a small value for \(m\), we can get a result similar to Laplace but with much lower cost.
posterior = calculate_posterior_VL(z,vecchia.approx,likelihood_model=data.model,
covparms = covparms)
if (spatial.dim==1){
par(mfrow=c(1,1))
ord = order(locs) # order so that lines appear correctly
y_limits = c(min(y, posterior$mean[ord]), max(y, posterior$mean[ord]))
plot(locs[ord], y[ord], type = "l", ylim = y_limits )
lines(locs[ord], posterior$mean[ord], type = "l", col=3, lwd=3)
legend("bottomright", legend = c("Latent", "VL"), col= c(1,3), lwd=c(1,3))
} else if (spatial.dim==2){
dfpar = par(mfrow=c(1,2))
# ordering unnecessary; we are using a scatter plot rather than lines
quilt.plot(locs, y, main= "Truth")
quilt.plot(locs, posterior$mean, main= "VL m=10")
par(defpar)
}
