  Profile likelihoods are rarely used in geostatistical models due to the
computational burden imposed by repeated decompositions of large variance
matrices. Accounting for uncertainty in covariance parameters can be highly
consequential in geostatistical models as some covariance parameters are poorly
identified, the problem is severe enough that the differentiability parameter
of the Matern correlation function is typically treated as fixed. The problem
is compounded with anisotropic spatial models as there are two additional
parameters to consider. In this paper, we make the following contributions: 1,
A methodology is created for profile likelihoods for Gaussian spatial models
with Mat\'ern family of correlation functions, including anisotropic models.
This methodology adopts a novel reparametrization for generation of
representative points, and uses GPUs for parallel profile likelihoods
computation in software implementation. 2, We show the profile likelihood of
the Mat\'ern shape parameter is often quite flat but still identifiable, it can
usually rule out very small values. 3, Simulation studies and applications on
real data examples show that profile-based confidence intervals of covariance
parameters and regression parameters have superior coverage to the traditional
standard Wald type confidence intervals.
