# Profile likelihoods for parameters in trans-Gaussian geostatistical models

Ruoyong Xu
Department of Statistical Sciences
University of Toronto
700 University Ave., Toronto, ON M5G 1Z5, Canada
<EMAIL_ADDRESS>
&Patrick Brown
Department of Statistical Sciences
University of Toronto
700 University Ave., Toronto, ON M5G 1Z5, Canada
<EMAIL_ADDRESS>

###### Abstract

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érn 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érn 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.

_Keywords_ profile likelihood $\cdot$ linear geostatistical models $\cdot$
parallel computing

##  1 Introduction

Anisotropic spatial models with the Matérn family of correlation functions are
important but difficult to implement because they involve many model
parameters. In addition to the five correlation parameters range ($\phi_{X}$),
shape ($\kappa$), nugget ($\nu^{2}$), anisotropy ratio ($\phi_{R}$),
anisotropy angle ($\phi_{A}$), there are the variance parameter ($\sigma$),
the Box-Cox transformation parameter ($\lambda$) and $p$ regression parameters
$\beta=(\beta_{1},\dots,\beta_{p}$) to estimate. Profile likelihood can be
used to reduce the dimensionality of the full likelihood, and it is generally
advised to plot the profile likelihood functions for intractable parameters.
In spatial models, the likelihoods of the correlation parameters can be far
from Normal (e.g., the profile likelihood surface of $\kappa$ can often be
very flat), and the true values of the shape, nugget or variance parameter can
sometimes be close to a boundary (e.g., true value of $\kappa$ can be
infinity). In these situations the curvature (related to the Fisher
information matrix) of the likelihood at the maximum likelihood estimates
(MLEs) is not a useful approximation of the rest of the likelihood surface,
and standard Wald confidence intervals (CIs) are inaccurate.

Most spatial software packages do not use profile likelihoods due to the
computational burden imposed by repeated likelihood evaluations. Let
$\omega=(\phi_{X},\kappa,\nu^{2},\phi_{Y},\phi_{A})$ be the vector of
correlation parameters, let $\psi=(\omega,\sigma,\lambda)$ and
$\beta_{i},i\in\\{1,p\\}$ denote an element of $\beta$. The geostatistics R
package _geostatsp_ fixes the parameters $\psi$ at their MLEs, denoted by
$\hat{\psi}$, and uses the pseudo-likelihood $\ell(\beta_{i};\hat{\psi})$ to
evaluate each $\beta_{i}$. For confidence intervals (CIs) of the covariance
parameters, second derivative of the profile likelihood
$\ell_{p}(\psi,\hat{\beta}_{\psi};y)$ is taken with respect to $\psi$
numerically to obtain the information matrix evaluated at $\hat{\psi}$.
Likelihood-based CIs would be advantageous in geostatistical context as the
‘boundary problem’ and ‘not regular likelihood’ often occurs. Obtaining
likelihood-based CIs is computed by maximizing over $\psi$, that is computing
the profile likelihood
$\ell_{p}(\beta_{i})=\ell(\beta_{i},\hat{\psi}_{\beta_{i}})$ and then deriving
CIs from the profile likelihood curves for each $\beta_{i}$. There is a
function in _geostatsp_ that is able to compute likelihood-based CIs for some
of the model parameters ($\lambda$ and $\omega$ only). For a fixed value of
$\phi_{X}$, for example, _geostatsp_ numerically optimize the joint likelihood
of the remaining parameters $(\kappa,\nu^{2},\phi_{R},\phi_{A})$. The process
is repeated for a large number of values of $\phi_{X}$, then it uses the
interpolated profile likelihood obtained by interpolating to compute
likelihood-based CIs for $\phi_{X}$. This procedure requires running an
optimization function a great many of times so can be very time-consuming, and
can fail if the spatial data has a large number of observations. Obtaining
plug-in Wald CIs is the relatively ‘easy way’, as quantities (MLEs and
information matrix) needed to be calculated are all available from the
optimizer used to find the MLEs of the correlation parameters. However,
$\ell(\beta_{i};\hat{\psi})$ does not allow for uncertainty in the estimated
correlation parameters to be reflected in the CIs for $\beta_{i}$’s, hence
_geostatsp_ ’s CIs for regression parameters would be theoretically narrower
than the likelihood-based CIs of $\beta_{i}$’s. Accounting for uncertainty in
correlation parameters is highly consequential in geostatistical models as
some correlation parameters are poorly identified, the problem is severe
enough that the differentiability parameter $\kappa$ of the Matérn correlation
is typically treated as fixed at some plausible value.

In this paper we present a method that is able to compute the profile
likelihoods and likelihood-based CIs for all the parameters in Gaussian linear
geostatistical models (LGMs) (anisotropic or isotropic), with dense covariance
matrices (Matérn), without approximating the covariance matrix. Our method
largely reduces the computation burden of profile likelihoods by taking
advantage of the parallel computing power of graphic processing units (GPUs),
so that the repetitive model fitting at a large number of different parameter
values can be done parallely on GPUs, and expedited by use of GPU’s local
memory. The method is implemented and made available in the R package
_gpuLik_.

The rest of this paper is organized as follows. In section 2, the necessary
background of our method is provided. Section 3 describes the main idea and
procedure and computational work of our method in detail. In Section 4 we
discuss briefly the considerations in GPU computation, and also explained how
the computational work is parallelized on a GPU. Section 5 presents 2
simulation studies. In section 6 we demonstrate the applicability of our
method by employing it on small and large real data sets. Section 7 concludes
with a discussion and future directions. Description on likelihood inference
using restricted maximum likelihood for linear geostatistical models is
deferred to the Appendix.

##  2 Preliminaries

The fit of the linear geostatistical model can often be improved by applying a
transformation to the observations $Y_{i},i=1,\dots,n$. Skewed random fields
that deviate from the Normal distribution can be modeled by means of the Box-
Cox transformation (BCT) [1]. For positive-valued observations $Y_{i}$, the
Box-Cox transformed response $Y_{i}^{\prime}$ has the form

$Y_{i}^{\prime}=b(Y_{i};\lambda)=\begin{cases}(Y_{i}^{\lambda}-1)/\lambda:&\lambda\in{\mathbb{R}}\quad\text{and}\quad\lambda\neq
0,\\\ \log(Y_{i}):&\lambda=0.\end{cases}$

The aim of the BCT is to make data closely resemble the Normal distribution,
so that the usual assumptions for linear models hold. The BCT cannot guarantee
that any probability distribution can be normalized. However, even in cases
when an exact mapping to the Normal distribution is not possible, BCT may
still provide a useful approximation [2]. Consider a zero-mean Gaussian random
field $U=\left[U(s_{1}),\dots,U(s_{n})\right]^{{{}^{{}_{\sf T}}}}$, the Box-
Cox transformed responses $Y^{\prime}_{i},i=1,\dots,n$ obtained at locations
$s_{i}$ with explanatory variables $X(s_{i})$ are modeled as

$\displaystyle
Y^{\prime}_{i}|U(s_{i})\overset{ind}{\sim}\text{N}(\mu(s_{i}),\tau^{2}),$
$\displaystyle\mu(s_{i})=X(s_{i})\beta+U(s_{i}).$

Given $U(s_{i})$, $Y_{i}^{\prime}$’s are mutually independently distributed
with Normal distribution with observation variance $\tau^{2}$.
$\beta=(\beta_{1},\dots,\beta_{p})^{{{}^{{}_{\sf T}}}}$ is the vector of
regression parameters. The covariance between locations $s+h$ and $s$ takes
the form below, where $\rho(\cdot)$ is the correlation function, and
$\sigma^{2}$ is the spatial variance or variability in residual variation.

$\displaystyle{\mathrm{COV}}\left[U(s+h),U(s)\right]=\sigma^{2}\rho(h;\phi,\kappa)$

Let $Y^{\prime}=(Y^{\prime}_{1},\dots,Y^{\prime}_{n})^{{{}^{{}_{\sf T}}}}$ be
the vector of transformed observations, and
$X=\left[X(s_{1}),\dots,X(s_{n})\right]^{{{}^{{}_{\sf T}}}}$ be an $n\times p$
design matrix, we rewrite the model in compact format

$\displaystyle
Y^{\prime}\sim\text{MVN}\left(X\beta,\sigma^{2}R(\phi,\kappa)+\tau^{2}I\right),$
$\displaystyle R_{ij}=\rho(s_{i}-s_{j};\phi,\kappa),$

here $R(\cdot)$ is the correlation matrix. Write $V=R(\phi,\kappa)+\nu^{2}I$,
and $\nu^{2}={\tau^{2}}/{\sigma^{2}}$, then

$Y^{\prime}\sim\text{MVN}\left(X\beta,\sigma^{2}V\right).$ (1)

###  2.1 Matérn correlation and Geometric anisotroy

Stein [3] strongly recommended using the Matérn class of spatial correlations
as this class has a parameter which controls the differentiability of the
process. We use the following parametrization of the Matérn, which includes a
rotation parameter to model directional effects.

$\displaystyle\rho(x;\phi,\kappa)=$
$\displaystyle\frac{2^{\kappa-1}}{\Gamma(\kappa)}\left[\sqrt{8\kappa}d(x)\right]^{\kappa}K_{\kappa}\left[\sqrt{8\kappa}d(x)\right]$
(2) $\displaystyle d(x)=$
$\displaystyle\left|\left|\begin{pmatrix}1/\phi_{X}&0\\\
0&1/\phi_{Y}\end{pmatrix}\cdot\begin{pmatrix}\cos\phi_{A}&-\sin\phi_{A}\\\
\sin\phi_{A}&\cos\phi_{A}\end{pmatrix}\cdot\begin{pmatrix}x_{1}\\\
x_{2}\end{pmatrix}\right|\right|$ (3)

Here $K_{\kappa}(\cdot)$ denotes the modified Bessel function of the second
kind of order $\kappa$, $\kappa$ is the shape parameter that determines the
differentiability of $U(x)$. For $\kappa=0.5$, the Matérn correlation function
coincides with the exponential correlation $\exp(-d)$. When
$\kappa\rightarrow\infty$, it tends to the Gaussian correlation
$\exp(-2d^{2})$. The standard isotropic correlation has a (scalar) range
parameter $\phi$ which controls how quickly correlation between two locations
falls as the distance between them increases. The anisotropic model has three
parameters which determine distance between points. The rotation parameter
$\phi_{A}$ gives the direction the correlation is strongest, and the two range
parameters $\phi_{X}$ and $\phi_{Y}$ which control the speed of decline in the
strongest and weakest direction. The isotropic model results when
$\phi_{X}=\phi_{Y}$, in which case the angle $\phi_{A}$ has no effect.

###  2.2 Inference based on profile likelihoods

Let $y^{\prime}=(y^{\prime}_{1},\dots,y^{\prime}_{n})^{{{}^{{}_{\sf T}}}}$ be
the Box-Cox transformation of the observations
$y=(y_{1},\dots,y_{n})^{{{}^{{}_{\sf T}}}}$. Further, write
$\omega=(\phi_{X},\kappa,\nu^{2},\phi_{Y},\phi_{A})$ as the vector of
covariance parameters. The log-likelihood function for the LGM parameters
given $y$ is

$\displaystyle-2\ell(\beta,\sigma^{2},\lambda,\omega;y)=$
$\displaystyle(y^{\prime}-X\beta)^{{{}^{{}_{\sf
T}}}}(\sigma^{2}V)^{-1}(y^{\prime}-X\beta)+\log|\sigma^{2}V|$
$\displaystyle-2*(\lambda-1)\sum_{i=1}^{n}\log y_{i}+n\log(2\pi),$ (4)

the third term on the right-hand side of (4) is the Jacobian of the Box-Cox
transformation.

Given $\omega$ and $\lambda$, The log-likelihood function (4) is maximized at

$\displaystyle\hat{\beta}_{\omega,\lambda}=(X^{{{}^{{}_{\sf
T}}}}V^{-1}X)^{-1}X^{{{}^{{}_{\sf T}}}}V^{-1}y^{\prime},\quad\text{and}$ (5)
$\displaystyle\hat{{\sigma}}^{2}_{\omega,\lambda}=\frac{1}{n}(y^{\prime}-X\hat{\beta}_{\omega,\lambda})^{{{}^{{}_{\sf
T}}}}V^{-1}(y^{\prime}-X\hat{\beta}_{\omega,\lambda}).$ (6)

Substituting (5) and (6) into (4), we obtain the profile log-likelihood (PLL)
for $\omega$ and $\lambda$

$\displaystyle-2\ell_{p}(\omega,\lambda;y)=$ $\displaystyle
n\log\frac{(y^{\prime}-X\hat{\beta}_{\omega,\lambda})^{{{}^{{}_{\sf
T}}}}V^{-1}(y^{\prime}-X\hat{\beta}_{\omega,\lambda})}{n}+\log|V|$
$\displaystyle-2(\lambda-1)\sum_{i=1}^{n}\log y_{i}+n\log(2\pi)+n.$ (7)

There is no closed formula for the MLE’s of the model parameters. Numerical
optimization of (2.2) yields the MLE’s $\hat{\omega}$ and $\hat{\lambda}$, and
back substituting them into (5) and (6) leads to the MLE’s $\hat{\beta}$ and
$\hat{\sigma}$. Inference for the model parameters using restricted maximum
profile likelihood is briefly presented in the Appendix.

###  2.3 Information-based and likelihood-based confidence intervals

We use $\hat{\theta}$ to denote the MLE of an arbitrary parameter vector
$\theta=(\theta_{1}\ldots\theta_{q})$. Assymptotic $100(1-\alpha)\%$
confidence intervals for a single parameter $\theta_{i}$ can be constructed in
two ways.

1\. Wald, or information-based confidence interval

$\\{\theta_{i}:|\theta_{i}-\hat{\theta_{i}}|\leq
z_{1-\alpha/2}[I(\hat{\theta})^{-1}]_{i,i}^{1/2}\\}$ (8)

Here $I(\hat{\theta})$ is the observd Fisher information evaluated at the MLE
and $z_{1-\alpha/2}$ is the $1-\alpha/2$ percentile of the standard Normal
distribution.

2\. Likelihood-based confidence interval

$\\{\theta_{i}:\ell[\theta_{i},\hat{\theta}_{-i}(\theta_{i})]-\ell(\hat{\theta})\geq
c_{1,\alpha}/2\\},\qquad\hat{\theta}_{-i}(\theta_{i})=\text{argmin}_{\theta_{-i}}\ell(\theta_{i},\theta_{-i})$
(9)

Here $\theta_{-i}$ refers to the vector $\theta$ without the $i$th element and
$(\theta_{i},\hat{\theta}_{-i})$ is referred to as the constrained MLE. The
$c_{1,\alpha}$ is a quantile of the Chi-squared distribution with
${\mathbb{P}}(\chi^{2}_{1}\leq c_{1,\alpha})=\alpha$.

Note that both results are asymptotic in $n$, the dimensionality of
$Y^{\prime}$. P. J. Diggle and Giorgi [4] and Haskard [5] are amongst those
who propose using likelihood-based CIs for spatial models, although both note
that the approximation does not hold if true value of $\theta$ is on a
boundary of the parameter space, or when one or more elements of $\theta$
define the range of $Y^{\prime}$. Baey, Courn‘ede, and Kuhn [6] and Visscher
[7] show that for non-spatial models when the variances of any subset of the
random effects are equal to zero, the asymptotic distribution of the test
statistics is a mixture of chi-square distributions with different degrees of
freedom. When true value of $\theta$ is away from the boundary the accuracy of
the approximation for geostatistical models has not been investigated.

##  3 Methodology

The main idea behind our method is to create likelihood-based confidence
intervals by evaluating the likelihood multiple times in parallel. We first
construct a quadratic approximation to the likelihood, and then to select a
large number of representative values for $\omega$ where we evaluate the
likelihood. An interpolant is built from these points and used to calculate
the profile likelihood. The time consuming part of the algorithm is evaluating
the likelihood many thousands of times, and our GPU implementation performs
computations for hundreds of parameter values simultaneously, with each
individual likelihood evaluation further parallelized and using thousands of
GPU work items.

The LGM has $7+p$ model parameters: five correlation parameters
$\omega=(\phi_{X},\kappa,\nu^{2},\phi_{Y},\phi_{A})$, the variance parameter
$\sigma^{2}$, Box-Cox parameter $\lambda$, and $p$ elements of $\beta$.
Profile likelihoods for all these parameters are computed using the following
steps.

  1. Step 1,

Obtain the MLE’s $\hat{\omega}$ and $\hat{\lambda}$ by numerical maximization
of the profile likelihood.

  2. Step 2,

Use a quadratic approximation of the log-likelihood to create a set of
representative parameters $\omega_{1},\dots,\omega_{K}$ and
$\lambda_{1},\dots,\lambda_{M}$.

  3. Step 3,

Evaluate and store various matrix determinants and cross products for the
$K\cdot M$ combinations of parameters $\omega_{k}$ and $\lambda_{k}$. Note
there are only $K$ different variance matrices.

  4. Step 4,

Compute profile likelihoods $\ell_{p}(\omega_{k},\lambda_{k})$ and obtain
univariate PLLs for correlation parameters and $\lambda$ using the highest
profile likelihoods.

  5. Step 5,

Obtain PLL for multiple values of $\sigma$ by taking the maximum of
$\ell_{p}(\sigma,\omega_{k},\lambda_{k},\hat{\beta}_{\sigma,\omega_{k},\lambda_{k}};y)$
over $k$ and $m$.

  6. Step 6,

Obtain PLL for $\beta$ by maximizing
$\ell_{p}(\beta,\omega_{k},\lambda_{k},\hat{\sigma}_{\beta,\omega_{k},\lambda_{k}};y)$
over $k$ and $m$. Section (3.4) shows how to get univariate PLL’s for
individual $\beta_{p}$ or pairs of coefficients.

###  3.1 Step 1: MLE’s of covariance parameters

Step 1 is carried out by numerically maximizing the profile likelihood
$\ell_{p}(\omega,\lambda;y)$ in (2.2). This requires sequential evaluations of
the likelihood and is done on the CPU using the L-BFGS-B algorithm implemented
in the _geostatsp_ package [8]. Since the likelihood is typically very flat
with respect to the shape parameter $\kappa$, two constrained maxima fixing
$\kappa=0.5$ and $\kappa=10$ are also computed.

###  3.2 Step 2: A set of representative parameters

A coarse configuration of the nuisance parameters will lead to bumps in the
profile curve [9], a large number of properly and sufficiently sampled
parameter values would be needed to make the profile curves smooth. In a
manner similar to the _R-INLA_ software package [10], we create a set of
“configuration" points $\omega_{k};k=1\ldots K$ using a quadratic
approximation to the profile likelihood $\ell_{p}(\omega,\lambda;y)$ in (2.2).
These points distributed evenly around the surface specified by
$\ell_{p}(\omega,\lambda;y)=c$, where $c$ takes on various $\chi_{2}$
quantiles from the likelihood’s maximum. If the approximation were exact, the
likelihood evaluated at $\omega_{k},\lambda_{k}$ would be equal to the
$\chi^{2}$ quantiles of the likelihood of the model parameters. We find a
reparameterization which makes exact and approximate likelihoods as similar as
possible.

Step 2 involves three components: reparametrizing the likelihood function;
computing a quadratic approximation using the Hessian matrix; and identifying
representative points.

_Reparameterization_

The model parameters are reparametrized in order to render the profile
likelihood as close to quadratic as possible, and thereby improving the
accuracy of the approximation. The transformation of the Matérn shape
parameter $\kappa$ depends on whether the MLE $\hat{\kappa}$ is large. The
reparametrization used is

$\tilde{\kappa}=\begin{cases}\log\kappa,&\hat{\kappa}<4\\\
1/\sqrt{\kappa},&\hat{\kappa}\geq 4.\end{cases}$

When $\hat{\kappa}$ is large the likelihood tends to be flat and the inverse-
root transform proved better than the log-transform in practice.

The three parameters controlling the rate of decay in the spatial correlation,
$\phi_{X}$, $\phi_{Y}$ and $\phi_{A}$ are reparametrized as proposed by [11]
with

$\displaystyle\gamma_{1}=\log\phi_{X}+\log\phi_{Y},\quad\text{(sumLogRange)}$
$\displaystyle\gamma_{2}=\sqrt{\phi_{X}/\phi_{Y}-1}\cos(2\phi_{A}),\quad\text{(aniso1)}$
$\displaystyle\gamma_{3}=\sqrt{\phi_{X}/\phi_{Y}-1}\sin(2\phi_{A}).\quad\text{(aniso2)}$

Note that $\exp(\gamma_{1}/2)$ (combinedRange) is the geometric mean of
$\phi_{X}$ and $\phi_{Y}$. The new parameterisation of
$(\gamma_{2},\gamma_{3})$ applies a polar coordinate transformation to the
ratio of ranges $\phi_{X}/\phi_{Y}$ and the angle $\phi_{A}$, which [11]
chooses because the isotropic model corresponds to the single point
$\gamma_{2}=\gamma_{3}=0$, avoiding the problem that the likelihood is flat in
the $\phi_{A}$ dimension when data are isotropic. The motivation for this
reparameterization can be more easily understood through results presented in
Section 6 and 7.

We call the parameters $\omega=(\phi_{X},\kappa,\nu^{2},\phi_{Y},\phi_{A})$
the ‘natural parameters’, and the reparameterized parameters
$\omega^{\prime}=(\gamma_{1},\tilde{\kappa},\nu,\gamma_{2},\gamma_{3})$ the
‘internal parameters’.

_Hessian matrix_

There are six non-linear parameters, which are the Box-Cox parameter $\lambda$
and the five covariance parameters in $\omega^{\prime}$. For creating a
quadratic approximation to the log-likelihood, we implicitly assume the second
derivatives of the log likelihood with respect to $\lambda$ and
$\omega^{\prime}$ are zero, and focus on the 5-dimensional Hessian matrix of
$\ell_{p}(\omega^{\prime},\hat{\lambda};y)$, varying $\omega^{\prime}$ but
fixing $\lambda$ at the MLE (as a function of $\omega^{\prime}$). The second
derivative of the likelihood with respect to $\lambda$ is used to approximate
the likelihood as a function of the Box-Cox parameter. Additionally, two
4-dimensional Hessians (fixing $\kappa$ at 0.5 and 10) are computed.

All Hessians are evaluated using numerical differentiation. For the
5-dimensional Hessian, we need to evaluate $51\cdot 3$ likelihoods
($\hat{\omega}^{\prime}\pm\delta$ and $\hat{\lambda}\pm\delta$), and for the
4-dimensional Hessian there are $33\cdot 3$ likelihoods to calculate. The
likelihoods are computed on GPU in parallel. Flat likelihood are fairly common
in spatial statistics,

[12] shows analytically that the range, shape and variance parameters are
often not well separately identifiable. It is possible that the Hessian is not
negative definite or with eigenvalues very close to zero. We take eigen-
decomposition of (minus) Hessian matrix $-H(\hat{\omega}^{\prime})$ (the
results will be required in obtaining representative points later):

$-H(\hat{\omega}^{\prime})=EDE^{{{}^{{}_{\sf T}}}},$

where $E$ is a matrix of eigenvectors and $D$ a diagonal matrix of
eigenvalues, and replace any negative eigenvalues with their absolute values.
If the range of the eigenvalues is large, i.e., the largest one is over 100,
while the smallest one is close to 0, then we replace the smallest eigenvalues
0.1, so to avoid one axes of the ellipsoid dominating the points selection.
This treatment of small and negative eigenvalues has worked well in practice
for a number of real and simulated datasets.

_Obtaining representative points_

Well-behaved log-likelihood surfaces can be well approximated around
$\hat{\omega}^{\prime}$ by a quadratic function derived from Taylor’s
expansion

$\ell(\omega^{\prime})=\ell(\hat{\omega}^{\prime})+\frac{1}{2}(\omega^{\prime}-\hat{\omega}^{\prime})^{{{}^{{}_{\sf
T}}}}(-H(\hat{\omega}^{\prime}))^{-1}(\omega^{\prime}-\hat{\omega}^{\prime}).$
(10)

The solid ellipsoid of $\omega^{\prime}$ values satisfying

$(\omega^{\prime}-\hat{\omega}^{\prime})^{{{}^{{}_{\sf
T}}}}(-H(\hat{\omega}^{\prime}))^{-1}(\omega^{\prime}-\hat{\omega}^{\prime})\leq\chi^{2}_{p}(\alpha)$
(11)

defines a confidence region with probability $1-\alpha$, where
$\chi^{2}_{p}(\alpha)$ is the upper $(100\alpha)$th percentile of a chi-square
distribution with $p$ degrees of freedom [13]. In order to get points that are
regularly-spaced on a hyper-ellipsoid surface, we first find points
$x=(x_{1},\dots x_{n})$ that are uniformly spaced on a 5-dimensional unit
hyper-sphere centered at the origin, by using a stochastic optimizer which
maximize the distance between the closest pair of points
$\max_{X}\min_{i,j}\|x_{i}-x_{j}\|$. We then transform the points $x$ to
$x^{\prime}$ on the desired hyper-ellipsoid surface by the formula

$\displaystyle x^{\prime}=E^{-2}D^{1/2}x+\hat{\omega}^{\prime}.$

There are 726 representative points ($n=726$) selected from each density
contour level. Parameter sets from usually 9 to 11 contour levels are combined
together to become the final configuration matrix for correlation parameters.
For $\lambda$, we set a vector of often 30 equally-spaced values between the
0.01 and 0.99 quantiles of its 1-dimensional quadratic approximation at
$\hat{\lambda}$.

Since the shape parameter is empirically difficult to estimate, a quadratic
approximation may not work well on the shape parameter when its log-likelihood
is far from quadratic. Hence, we fit several additional models with the Matérn
shape parameter fixed at values between 0.5 and 20. The Hessian matrix is for
these additional models is 4 by 4, with $n=120$ representative points from
each density contour level. The likelihood is often very flat for the shape
parameter, this procedure let us explore the likelihood for a range of
$\kappa$ values. It is possible to sample negative values for nugget $\nu^{2}$
from the quadratic approximation if the true nugget is close to 0, in that
case, we turn half of the negative values to 0, and the other half to random
uniform numbers sampled from the interval $(0,2)$, see Figure 7(c) for a plot
of the sampled nugget values and their corresponding profile log-likelihoods
in a real data example, in which a visible small pile of points is collected
at 0.

###  3.3 Step 3: Log-likelihood computation

The computationally expensive steps in evaluating the likelihood are taking
the Cholesky decomposition of the variance matrix $V$ and performing
triangular backsolve computations. These steps are performed on the GPU, and
the crossproducts and matrix determinants listed in Table 1 are saved and
returned to the CPU. These components of the likelihoods are then used to
compute different variations of the likelihood (i.e. profiling out $\beta$ or
$\sigma$ or both).

Component name in R $\log|V|$ detVar $\log|X^{{{}^{{}_{\sf T}}}}V^{-1}X|$
detReml $y^{{{}^{{}_{\sf T}}}}V^{-1}y$ ssqYX upper-left diagonal entries
$X^{{{}^{{}_{\sf T}}}}V^{-1}X$ ssqYX lower-right diagonal entries
$X^{{{}^{{}_{\sf T}}}}V^{-1}y$ ssqYX lower-left block
$(X\hat{\beta})^{{{}^{{}_{\sf T}}}}V^{-1}(X\hat{\beta})$ ssqBetahat
$(y-X\hat{\beta})^{{{}^{{}_{\sf T}}}}V^{-1}(y-X\hat{\beta})$ ssqResidual

Table 1: Summary statistics computed and retained for each set of parameter
values

We take the following steps to compute the items in Table 1 for the parameter
sets $\omega^{\prime}_{1}$ to $\omega^{\prime}_{K}$.

  1. Step 1:

Compute the Matérn correlation matrices $V=(V_{1},\dots,V_{K})$ with the
maternBatch() function.

  2. Step 2:

Take the Cholesky decompositions $V_{k}=L_{k}D_{k}L_{k}^{{{}^{{}_{\sf T}}}}$
with cholBatch(), where the $L_{k}$ are lower unit triangular matrices and the
$D_{k}$ are the diagonal matrices. The function, then
$\log|V|=2\log|L|+\log|D|=\log|D|$. cholBatch() returns $L_{k}$, $D_{k}$ and
$\log|V_{k}|=\text{sum}(\log D_{k})$.

  3. Step 3:

Solve $L_{k}^{-1}(y,X)$ with backsolveBatch().

  4. Step 4:

Use the above to compute $(y,X)^{{{}^{{}_{\sf T}}}}V_{k}^{-1}(y,X)$ with
crossprodBatch(), the results are saved in ssqYX.

  5. Step 5:

Perform a Cholesky decomposition on $X_{k}^{{{}^{{}_{\sf
T}}}}V_{k}^{-1}X_{k}=Q_{k}P_{k}Q_{k}^{{{}^{{}_{\sf T}}}}$ with cholBatch() and
calculate $\log|X_{k}^{{{}^{{}_{\sf T}}}}V_{k}^{-1}X_{k}|$.

  6. Step 6:

Solve the equation $c=Q_{k}^{-1}X^{{{}^{{}_{\sf T}}}}V_{k}^{-1}y$ with
backsolveBatch().

  7. Step 7:

Compute $\text{\tt ssqBetahat}=(X\hat{\beta})^{{{}^{{}_{\sf
T}}}}V^{-1}(X\hat{\beta})=c^{{{}^{{}_{\sf T}}}}P_{k}^{-1}c$ by
crossprodBatch().

  8. Step 8:

Compute $\text{\tt ssqResidual}=(y-X\hat{\beta})^{{{}^{{}_{\sf
T}}}}V^{-1}(y-X\hat{\beta})=\text{\tt ssqYX}-\text{\tt ssqBetahat}$.

GPU memory limitations mean not all $K$ of the parameter sets can be processed
simultaneously, the limiting factor is the memory occupied by the variance
matrices $V_{k}$. The results presented here process several hundred parameter
sets at a time, the maximum number possible depends on the size of the dataset
and the amount of GPU memory.

###  3.4 Profile likelihood for $\beta_{p}$

Here we consider computing the profile likelihood for a single regression
coefficient $\beta_{p}$ for the $p$th covariate. The profile likelihood is

$\displaystyle\ell_{p}(\beta_{p};y)$
$\displaystyle=\underset{\beta_{-p},\sigma,\omega,\lambda}{\text{max}}\\{\ell(\beta_{p},\beta_{-p},\sigma,\omega,\lambda;y)\\}$
$\displaystyle\approx\underset{k,m}{\text{max}}\\{\ell(\beta_{p},\hat{\beta}_{-p}(\beta_{p},\omega_{k},\lambda_{\ell}),\hat{\sigma}(\beta_{p},\omega_{k},\lambda_{\ell}),\omega_{k},\lambda_{m})\\}$
(12)

where $\beta_{-p}$ is the vector $\beta$ without its $p$th element. The
functions $\hat{\sigma}(\cdot)$ and $\hat{\beta}_{-p}(\cdot)$ denote the the
MLE’s of $\beta_{-p}$ and $\sigma$ for given values of $\beta_{p}$, $\omega$
and $\lambda$, for which there are closed-form solutions.

When $\beta_{p}$ is fixed, the term $X_{p}\beta_{p}$ becomes a fixed offset
and the term in the model involving covariates becomes $X_{-p}{{{}^{{}_{\sf
T}}}}\beta_{-p}$, where $X_{-p}$ is the matrix $X$ without the $p$th column.
Computing $\hat{\beta}_{-p}(\beta_{p},\omega,\lambda)$ involves the following
matrix crossproducts:

$\displaystyle(y-X_{p}\beta_{p})^{{{}^{{}_{\sf T}}}}V^{-1}(y-X_{p}\beta_{p})$
$\displaystyle=y^{{{}^{{}_{\sf T}}}}V^{-1}y-2\beta_{p}(X^{{{}^{{}_{\sf
T}}}}V^{-1}y)_{[p,]}+\beta_{p}^{2}(X^{{{}^{{}_{\sf T}}}}V^{-1}X)_{[p,p]},$
$\displaystyle(X_{-p}\hat{\beta_{-p}})^{{{}^{{}_{\sf
T}}}}V^{-1}(X_{-p}\hat{\beta_{-p}})$ $\displaystyle=(X^{{{}^{{}_{\sf
T}}}}V^{-1}y)^{{{}^{{}_{\sf T}}}}_{[-p,]}(X^{{{}^{{}_{\sf
T}}}}V^{-1}X)^{-1}_{[-p,-p]}(X^{{{}^{{}_{\sf T}}}}V^{-1}y)_{[-p,]}$
$\displaystyle-2\beta_{p}(X^{{{}^{{}_{\sf T}}}}V^{-1}X)^{{{}^{{}_{\sf
T}}}}_{[-p,p]}(X^{{{}^{{}_{\sf T}}}}V^{-1}X)^{-1}_{[-p,-p]}(X^{{{}^{{}_{\sf
T}}}}V^{-1}y)_{[-p,]}$ $\displaystyle+\beta_{p}^{2}(X^{{{}^{{}_{\sf
T}}}}V^{-1}X)^{{{}^{{}_{\sf T}}}}_{[-p,p]}(X^{{{}^{{}_{\sf
T}}}}V^{-1}X)^{-1}_{[-p,-p]}(X^{{{}^{{}_{\sf T}}}}V^{-1}X)_{[-p,p]}.$

The subscript notation $[-p,-p]$ indicates the $p$th row and column is
excluded from a matrix, with $[-p,p]$ refers to the $p$th column of a matrix
while excluding the $p$th row. Computing the above only needs subsetting rows
and columns from the ssqYX matrices listed in Table 1, and the profile
likelihoods for each $\beta_{p}$ can be computed without additional use of the
GPU.

###  3.5 Profile likelihood for $\sigma$

By substituting (5) into (4), it is easy to obtain the PLL for
$(\sigma,\omega,\lambda)$ as

$\displaystyle-2\ell_{p}(\sigma,\omega,\lambda;y)-n\log(2\pi)=$
$\displaystyle(y-X\hat{\beta})^{{{}^{{}_{\sf
T}}}}(\sigma^{2}V)^{-1}(y-X\hat{\beta})+n\log\sigma^{2}$
$\displaystyle+\log|V|-2(\lambda-1)\sum_{i=1}^{n}\log y_{i},$ (13)

where $\hat{\beta}$ and $V$ both depend on $\omega$. The various terms in
(3.5) can all be computed with the quantities in Table 1 and the exception of
the final term (the Box-Cox Jacobian) which is also returned by the GPU. We
approximate the PLL for $\sigma$ with

$\ell_{p}(\sigma)\approx\max_{k,m}\ell_{p}(\sigma,\omega_{k},\lambda_{m};y),$

and evaluate $\ell_{p}(\sigma)$ for a range of values of $\sigma$.

###  3.6 Profile likelihood for correlation parameters

We obtain the profile likelihood for any of the correlation parameters, either
reparametrized
$\omega^{\prime}=(\gamma_{1},\tilde{\kappa},\nu,\gamma_{2},\gamma_{3})$ or
original $\omega=(\phi_{X},\kappa,\nu^{2},\phi_{Y},\phi_{A})$, from the points
cloud of all computed likelihoods $\ell_{p}(\omega_{k},\lambda_{m};y)$ like
the example plots shown in Figure 7. To compute the PLL for the shape
parameter $\kappa$, for example, we create a 2D scatterplot from the
configuration points $\omega_{k}$ with $\kappa_{k}$ on the x-axis and
$\underset{m}{\text{max}}\ell_{p}(\omega_{k},\lambda_{m};y)$ on the y-axis.
The collection of ‘best’ likelihoods are obtained by computing the convex hull
of these points, shown as the 12 black circular points in Figure 7(b). The
green line in Figure 7(b), a linear interpolation of these 12 points, is the
estimated PLL curve. The PLL for the Box-Cox parameter $\lambda$ is simply the
interpolation of $\lambda_{m}$ and
$\underset{k}{\text{max}}\ell_{p}(\omega_{k},\lambda_{m};y)$.

Similarly, for 2-dimensional profile likelihood of any 2 of the correlation
parameters $(\omega^{\prime}_{i},\omega^{\prime}_{j})$ in $\omega^{\prime}$,
we first find the smallest convex hull that encloses the set of points
$[\omega^{\prime}_{i},\omega^{\prime}_{j},\underset{\lambda}{\text{max}}\\{\ell_{p}(\omega,\lambda;y)\\}]$
in the 3-dimensional space, remove points on the bottom facet of the hull, and
then the profile likelihood surface is linearly interpolated. We create
regular girds on the space spanned by
$(\omega^{\prime}_{i},\omega^{\prime}_{j})$, and predict values for these
grids using the interpolation.

For plotting 2-dimensional profile likelihoods of the correlation parameters
in the “natural space" $(\omega_{i},\omega_{j})$, we convert the regular girds
created in “natural space" to grids in “internal space", and prediction for
those grid points still uses interpolation formed in the “internal space". The
convex hull and interpolation method works poorly in the original parameter
space, since the log likelihood is less bell-shaped than in the reparametrized
space.

##  4 Computational considerations

Taking full advantage of the parallel computing abilities of GPU’s requires
some manual configuration of the GPU memory and work items (the GPU equivalent
of CPU cores). Work items are divided into work groups which share local
memory (which is faster than global memory). We treat work items as being on a
two-dimensional grid (a third dimension is unused), and divide tasks amongst
work items as shown in Table 2.

Work items are arrayed on an Nglobal0 by Nglobal1 grid, with each work group
containing Nlocal0 by Nlocal1 work items. Suppose Nglobal0=6, Nglobal1=4,
Nlocal0=3 and Nlocal1=2, which means we have $2\times 2$ work-groups on the
2-dimensional (dimension 0 and dimension 1) computation domain (See Figure 1),
each work-group consists of $3\times 2$ work-items. Suppose we have a total of
500 parameter configurations divided into 5 sets of 100. With each looping
inside the C++ function likfitGpu_Backend(), the GPU computes 100 log-
likelihoods and saves corresponding 100 sufficient statistics on the GPU’s
global memory.

Table 2: List of the main GPU kernels showing the tasks they compute and how
the task is divided amongst work items and work groups.

$\begin{array}[]{ lccccc
}\hline\cr\hline\cr\text{Kernel}&\text{Formula}&\text{group0}&\text{group1}&\text{local0}&\text{local1}\\\
\hline\cr\text{maternBatch}&V&\text{shared coordinates}&\text{shared
parameters}&\text{matrix element}&\text{matrix}\\\
\text{cholBatch}&V=LDL^{{{}^{{}_{\sf T}}}}&\text{matrix}&unused&\text{row of
L}&\text{column of L}\\\ \text{backsolveBatch}&AC=B&\text{matrix}&\text{col of
B and C}&\text{row of A and C}&\text{col of A, row of C}\\\
\text{crossprodBatch}&C=A^{{{}^{{}_{\sf T}}}}DA&\text{column of
C}&\text{matrix}&\text{elements of D}&\text{row of C}\\\
\hline\cr\hline\cr\end{array}$

Table 2 shows the four main GPU kernels which are used to perform matrix
operations. In each of the loop in this example, the _maternBatch_ kernel has
work items in dimension 0 looping through elements of the variance matrix and
work items in dimension 1 looping through parameter sets (each one used for a
different variance matrix). In the scenario above each work group will use
half of the 100 parameter sets and copy 50 sets of parameters to local memory,
since there are two work groups in dimension 1. Each work-group copies
NparamPerIter$/2=50$ parameters and 3 (the first element of Nlocal)
coordinates in the local memory. At each iteration a work group will copy over
pairs of coordinates for three matrix elements to local memory (since
Nlocal0=3), and use these coordinates for two matrices (since Nlocal1 = 2).

The Cholesky decomposition is performed by _cholBatch_ and only uses a single
work group in dimension 1 (requiring Nglobal1 = Nlocal1). The first work-group
takes Cholesky decomposition for variance matrices $V_{0},V_{2},V_{4},\dots$,
the second work-group takes cholesky decomposition for
$V_{1},V_{3},V_{5},\dots$. Within each work-group, the first-row work-items
are responsible for the $1,4,7,\dots$ rows of the first column of $L$, and
simultaneously, the second-row work-items is responsible for the $2,5,8,\dots$
rows of the first column of $L$, the third-row calculates cells in the
$3,6,9,\dots$ rows of the first column of $L$.

In the _backsolveBatch_ kernel, work-groups of the same row computes the same
$C_{i}$’s, in the example here, the upper two work-groups compute
$C_{0},C_{2},C_{4},\dots$ in order, and the lower two work-groups computes the
rest $C_{1},C_{3},C_{5},\dots$ in sequence. Each work-group is responsible for
ceil(ncol(B)/2) columns of B, in other words, work-groups shaded grey color
here are responsible for the columns $0,2,4,\dots$ of each $C_{i}$, and the
white work-groups are responsible for columns $1,3,5,\dots$ of each $C_{i}$.
Within each work-group, work-items on row 0 computes each of the rows
$0,3,6,\dots$ of each $C_{i}$ together, similarly and simultaneously, work-
items on row 1 computes each of the rows $1,4,7,\dots$ of each $C_{i}$
together, and work-items on row 2 computes each of the rows $2,5,8,\dots$ of
each $C_{i}$ together.

The _crossprodBatch_ kernel computes $C=A^{{{}^{{}_{\sf T}}}}DA$ or
$C=A^{{{}^{{}_{\sf T}}}}D^{-1}A$ or $C=A^{{{}^{{}_{\sf T}}}}A$. The two shaded
work-groups on the left are responsible for computing
$C_{O},C_{2},C_{4},\dots$, and the white work-groups are responsible for
computing $C_{1},C_{3},C_{5},\dots$. The group index in dimension 0 decides
the columns of C that the work-items in this work-group compute for, here the
upper work-group computes for the columns $0,2,4,\dots$ in each $C_{i}$, and
the lower work-group computes the columns $1,3,5,\dots$ of each $C_{i}$.
Within each work-group, work-items of same local index in dimension 1 compute
the same rows of $C_{i}$, so here the left three work-items compute rows
$0,2,4,\dots$ of each $C_{i}$ and the right three work-items compute rows
$1,3,5,\dots$ of each $C_{i}$.

Figure 1: Representation of 2-dimensional computation domain for the simple
example for illustration.

There is the argument NparamPerIter in the function likfitLgmGpu() which
specifies how many parameter sets are simultaneously computed on each
iteration of calls to the back-end GPU functions. The batch size value
(NparamPerIter) is limited by the available GPU global memory. If the LGM to
estimate has $n$ observations, $p$ covariates, and the vector configuration
for the Box-Cox parameter has $m$ values. Then each loop inside the backend
C++ function refreshes the GPU matrix $V$ (stores Matérn correlation
matrices), which has $n$*NparamPerIter rows and $n$ columns; _crossprodBatch_
’s result ssqYX is a GPU matrix of $(m+p)$*NparamPerIter rows and $m+p$
columns. We should make NparamPerIter relatively small if it is a large data
set to fit or if there are lots of covariates in the model and configuration
for $\lambda$ has a large number of values. The values of Nglobal and Nlocal
are better to be powers of two, and should not exceed the maximum number of
work-items on the GPU, which can be found by gpuR::gpuInfo().

##  5 Two simulation studies

To assess the CIs produced by our method in a realistic setting, we performed
2 simulation studies, study A uses anisotropic data and study B uses isotropic
data. Within each study, we do 150 simulations using simulated coordinates,
and the other 150 simulations using coordinates from the classic _Loaloa_
dataset.

Specifically, in study A, we generate 150 realizations $y_{1},\dots,y_{150}$
twice from the model $Y_{i}\sim N(5+1*X_{i1}+1*X_{i2}+U(s_{i}),\tau^{2})$,
where $\tau=0.8$. In the two simulations A(i) and A(ii), the stationary
Gaussian random field $U(s)$ with variance $\sigma^{2}=1$ and Matérn
correlation function with parameters $\kappa=2,\phi_{R}=2,\phi_{A}=0.2$ is
produced over different coordinates, and the two explanatory variables are
defined in different ways accordingly. No Box-Cox transformation was performed
on $Y$ so the parameter $\lambda=1$.

For study A(i), $U(x)$ is generated over 224 (Ngrid = 15) point locations on a
$9km\times 9km$ square (displayed in Figure LABEL:fig:sim1coords-1), the point
locations in Figure LABEL:fig:sim1coords-1 are randomly generated, we removed
one point that was very close to another point. The range parameter $\phi_{X}$
is set to be 1,000 for this coordinate. Covariates $X_{i1}$ and $X_{i2}$ were
defined using the coordinate values: $X_{i1}=coordsX_{i}/10000$,
$X_{i2}=(coordsX_{i}/10000)^{2}$, where $coordsX_{i}$ represents the x
coordinate of location $s_{i}$.

For study A(ii), $U(s)$ is generated on the _Loaloa_ coordinates, shown in
Figure LABEL:fig:sim1coords-2. _Loaloa_ is a real data set (available in the
_geostatsp_ package) that has 190 locations (from 190 village surveys), many
of them are clustered together. For A(ii), $\phi_{X}$ is set as 50,000.
Correspondingly, the two explanatory variables were defined as:
$X_{i1}=\log(elevationLoa_{i});X_{i2}=eviLoa_{i}*10^{-7}$, where
$elevationLoa_{i}$ means the log of elevation at location $s_{i}$.

Similarly, for study B, we produce a total of 300 realizations
$y_{1},\dots,y_{300}$ from the model $Y_{i}\sim
N(2+1*X_{i1}+1*X_{i2}+U(s_{i}),\tau^{2})$ with $\tau=0.8$ and $U(s)$ defined
on the same two types of coordinates as in study A. $U(s)$ have the same
parameter setup as in simulation A, except for $U(s)$ being isotropic. Again,
study B(i) has realizations produced on the simulated coordinates, study B(ii)
uses the _Loaloa_ coordinates.

(a)

(b)

Figure 2: Simulation coordinates for $U(s)$.

We computed coverage for 80% CIs for each of the 8 parameters and compared the
CI coverage rates for these parameters with results computed using _geostatsp_
in Table 3. The 150 PLL curves using simulated coordinates for the 8
parameters in study A(i) and B(i) are shown in Figure 3 and 4 respectively.
2-dimensional plots of the 80% confidence regions of the anisotropic
parameters are shown in Figure 5.

The results demonstrate that _gpuLik_ (or profile-likelihood-based CIs) have
higher and closer to expected coverage than _geostatsp_ ’s (or information-
based CIs). The superiority is much more obvious for the regression parameter
$\beta$’s and for simulation based on the real coordinates _Loaloa_.

For simulations using isotropic data, we also counted the coverage of the
2-dimensional 80% profile likelihood-based confidence regions of
$(\gamma_{2},\gamma_{3})$, which are shown in Figure 5. Out of 150
simulations, 103 confidence regions covered the true values of the parameters
$(\gamma_{2},\gamma_{3})$ (rate $=0.687$).

(a) $\log(\phi_{X}\phi_{Y})$ (b) $\kappa$ (c) $\nu^{2}$ (d) $\gamma_{2}$ (e)
$\gamma_{3}$ (f) Intercept (g) $\beta_{1}$ (h) $\beta_{2}$

Figure 3: Profile log-likelihoods plots for the model parameters in simulation
Study A (i) (Anisotropic data and simulated coordinates).

(a) $\log(\phi_{X}\phi_{Y})$ (b) $\kappa$ (c) $\nu^{2}$ (d) aniso1 (e) aniso2
(f) Intercept (g) $\beta_{1}$ (h) $\beta_{2}$

Figure 4: Profile log-likelihoods for the model parameters in simulation Study
B (i) (Isotropic data and simulated coordinates).

(a) $\gamma_{2}$ vs $\gamma_{3}$ (b) $\gamma_{2}$ vs $\gamma_{3}$

Figure 5: Plot of the 2-dimensional 80% confidence regions of the profile log-
likelihoods for $\gamma_{2}$ (on the x axis) and $\gamma_{3}$ (on the y axis)
from the 150 simulations in Study A (i) (left panel) using anisotropic data
and regular coordinates and B (i) (right panel) using isotropic data and
regular coordinates, respectively.  Table 3: A comparison between the coverage
rates of the 150 80% confidence intervals of the model parameters in
simulation studies A and B using methods provided by _gpuLik_ and _geostatsp_
respectively.

Anisotropic Isotropic Simulated coords Loaloa coords Simulated coords Loaloa
coords Parameters likelihood-based Wald likelihood-based Wald likelihood-based
Wald likelihood-based Wald Intercept 0.600 0.340 0.740 0.513 0.666 0.340 0.686
0.533 $\beta_{1}$ 0.800 0.786 0.813 0.173 0.800 0.786 0.813 0.180 $\beta_{2}$
0.800 0.766 0.820 0.180 0.793 0.746 0.813 0.266 $\phi_{X}$ or
$\sqrt{\phi_{X}\phi_{Y}}$ 0.760 0.460 0.760 0.706 0.420 0.193 0.513 0.400
$\kappa$ 0.967 0.800 0.926 0.766 0.980 0.393 0.933 0.506 $\nu^{2}$ or $\tau$
0.860 0.780 0.826 0.713 0.753 0.320 0.826 0.466 $\gamma_{2}$ or $\phi_{R}$
0.780 0.713 0.660 0.580 0.653 0.186 0.660 0.320 $\gamma_{3}$ or $\phi_{A}$
0.713 0.680 0.686 0.573 0.680 0.213 0.733 0.353 $\lambda$ 0.806 0.720 0.793
0.686 0.813 0.326 0.820 0.480

##  6 Example: Swiss rainfall data

Figure 6: Swiss rainfall data (colored blue points and top legend) with
elevation (background colors and bottom legend).

We now implement our method on the Swiss rainfall data from the _geostatsp_
package (See [8, 14] for a detailed description of the data and the SIC-97
project). The data, plotted in Figure 6, consists of 100 daily rainfall
measurements made in Switzerland in May, 1986. We fit the data with the Linear
Geostatistical Model

$\displaystyle Y^{\prime}_{i}\sim\text{N}[X(s_{i})^{{{}^{{}_{\sf
T}}}}\beta+U(s_{i}),\tau^{2}]$
$\displaystyle\text{COV}[U(s+h),U(s)]=\sigma^{2}\rho(h;\phi_{X},\phi_{Y},\phi_{A},\kappa)$

where $Y^{\prime}_{i}$ is the (Box-Cox transformed) rainfall measurement at
the $i$th location, the covariates $X(s)$ are comprised of an intercept and
the elevation at location $s$, and $U(s)$ is a zero-mean Gaussian process. The
correlation function $\rho$ is an anisotropic Matérn with range parameters
$\phi_{X}$ and $\phi_{Y}$, a rotation parameter $\phi_{A}$, and shape
parameter $\kappa$. Below the LGM is fit using swissRain and swissAltitude
objects from _geostatsp_ , with the geostatsp::lgm() function computing the
MLE’s of the model parameters.

    R> data("swissRain", package = "geostatsp")
    R> swissRain$elevation = raster::extract(swissAltitude, swissRain)
    R> mle1 = geostatsp::lgm(formula = rain ~ elevation, data = swissRain,
    +      covariates = swissAltitude, fixBoxcox = FALSE, fixShape = FALSE,
    +      fixNugget = FALSE, aniso = TRUE, reml = FALSE, grid = 20)

As described in Step 1 of the algorithm in Section 3, we fit five additional
models fixing the Matérn shape parameter $\kappa$ at 0.5, 0.9, 10, 20, and 100
(using arguments fixNuggeet = FALSE and shape = 0.5). A list containing these
six outputs from lgm() is passed to the gpuLik::configParams() function to
create representative points.

    R> reParamaters = gpuLik::configParams(list(mle1, mle2, mle3, mle4, mle5,
    +      mle6), alpha = c(0.00001, 0.01, 0.1, 0.2, 0.25, 0.3, 0.5, 0.8,
    +      0.9, 0.95, 0.99, 0.999))
    R> paramsUse <- reParamaters$representativeParamaters
    R> b <- reParamaters$boxcox

We select 12 significance levels for $\alpha$ between $0.00001$ and $0.999$
for the contours of the approximate likelihood surface, giving 15318 rows
representative points (including 6 MLEs from the 6 fitted LGMs).

Likelihoods for these representative parameters are evaluated by the
gpuLik::likfitLgmGpu() function below. The values for boxcox are a sequence of
33 equally-spaced values between the 0.01th and 0.99th quantiles of the Wald
CI, as well as the MLE $\hat{\lambda}$. The NparamPerIter argument specifies
that 400 likelihoods will be evaluated simultaneously, meaning a 40,000 by 100
matrix is needed to hold the 400 covariance matrices. The matrix
$\text{ssqYX}=(y,X)^{{{}^{{}_{\sf T}}}}V^{-1}(y,X)$ returned by
crossprodBatch() is of $36\cdot 400$ rows and $36$ columns,34 columns for the
34 Box Cox parameters and one column for each of the two covariates.

    R> result <- gpuLik::likfitLgmGpu(model = mle1, params = paramsUse,
    +    paramToEstimate =
    +      c(’range’, ’combinedRange’, ’shape’, ’nugget’, ’aniso1’, ’aniso2’, ’boxcox’),
    +  Ψboxcox = seq(b[1], b[9], len = 33), cilevel = 0.9,
    +    NparamPerIter = 400, Nglobal = c(128, 128), Nlocal = c(16, 16), NlocalCache = 2800)

(a) $\sqrt{\phi_{X}\phi_{Y}}$ or $\gamma_{1}$ (b) $\kappa$ (c) $\nu^{2}$ (d)
aniso1, $\gamma_{2}$ (e) aniso2, $\gamma_{3}$ (f) anisotropy ratio, $\phi_{R}$
(g) anisotropy radians, $\phi_{A}$ (h) Legend

Figure 7: Profile log-likelihoods for the correlation parameters estimated
from the Swiss rain data. Vertical axes are log likelihood minus the
likelihood at the MLE. Dashed lines show the MLE and 90% confidence intervals.
Point clouds are representative points and their likelihoods, with colours
showing the $\alpha$ quantile of the contour each point was drawn from.

The NlocalCache argument is related to the amount of local memory available on
the GPU, and 2,800 real numbers will be stored in local memory at any one
time. Each work group consists of a grid of 16 by 16 work items, comprising an
8 by 4 grid of work groups (giving $16\cdot 8$ by $16\cdot 4$ work items in
total).

Figure 7 contains profile likelihoods for covariance parameters. The point
clouds are the representative parameters and their log likelihoods (minus the
log likelihood at the MLE), with colours corresponding to the significance
level $\alpha$ of the contours the points were drawn from. The PLL curve for
$\kappa$ displayed in Figure 7(b) is quite flat on the right side of its
maximum, and the quadratic approximation to the likelihood is poor. The
configuration points with $\kappa$ fixed at at 0.5, 0.9, 10, 20, 100 are
clearly visible as vertical lines. The PLL for $\nu^{2}$ has a visible pile of
points collected at 0 as we turn half of the negative $\nu^{2}$’s in the
sample to be 0 and the other half are randomly sampled from the interval
$(0,2)$. The log-likelihoods for the reparametrized anisotropy parameters
$\gamma_{2}$ and $\gamma_{3}$ are fairly close to quadratic, which supports
the use of the polar coordinate transformation. The data are clearly
anisotropic although the ratio parameter has considerable uncertainty,
possibly being as small as 4 or as large as 15.

(a) Box-Cox, $\lambda$ (b) Standard deviation, $\sigma$ (c) Intercept,
$\beta_{0}$ (d) Elevation, $\beta_{1}\cdot 1000$

Figure 8: Profile log-likelihoods for the Box Cox transformation parameter,
standard deviation and regression parameters computed from the Swiss rain
data. Horizontal dashed lines on each plot define approximate 90% confidence
intervals.

The 2-dimensional PLL contour plots in Figure 9… of $(\gamma_{2},\gamma_{3})$
and $(\phi_{R},\phi_{A})$ in Figure 9 show the data are anisotropic. They
demonstrate the advantages of the reparameterization, the contours of the
internal parameters has a more oval-like shape compared to the PLL contours of
the natural parameters. As the value of anisotropy ratio gets closer to 1
(isotropic), the estimate of CI for the anisotropy radians would become wider
and eventually on the whole interval of $(-\pi/2,\pi/2)$.

(a) $(\beta_{0},\beta_{1}*1000)$ (b) $(\phi_{R},\phi_{A})$ (c)
$(\nu^{2},\gamma_{2})$ (d) $(\gamma_{1},\gamma_{2})$ (e)
$(\gamma_{2},\gamma_{3})$ (f) $(\gamma_{1},\gamma_{3})$

Figure 9: Two-dimensional profile log-likelihood contour plots for pairs of
model parameters, computed from the Swiss rain data. Circular plotting symbols
show locations of representative points on convex hull.

Table 4 shows the profile-based CIs from _gpuLik_ are wider than Wald based
CIs from _geostatsp_. The CIs for $\lambda$ obtained by the two methods are
very close.

Likelihood-based confidence intervals for model parameters are compared with
_geostatsp_ ’s Wald-based CI’s in Table 4. PLL plots for the model parameters
and some 2-dimensional PLL plots are displayed in Figures 7, 8 and 9.

Table 4: Swiss rainfall data. Comparison between Wald and profile likelihood based CI’s. | likelihood-based | Wald
---|---|---
| notation | estimate | ci0.05 | ci0.95 | ci0.05 | ci0.95
(Intercept) | $\beta_{1}$ | 5.01 | 2.76 | 7.8 | 3.32 | 6.69
elevation*1000 | $\beta_{2}*1000$ | 0.23 | -0.44 | 0.91 | -0.36 | 0.82
sdSpatial | $\sigma$ | 2.68 | 1.63 | 4.74 | 1.57 | 4.58
range/1000 | $\phi_{X}/1000$ | 38.62 | 21.79 | 118.63 | 18.54 | 80.42
combinedRange/1000 | $\sqrt{\phi_{X}\phi_{Y}}/1000$ | 13.58 | 6.34 | 48.96 | NA | NA
anisoRatio | $\phi_{R}$ | 8.09 | 3.68 | 14.58 | 4.73 | 13.86
shape | $\kappa$ | 1.83 | 0.40 | $>100$ | -1.25 | 4.91
nugget | $\nu^{2}$ | 0.14 | 0.00 | 0.28 | NA | NA
sdNugget | $\tau$ | 0.99 | 0.00 | 1.41 | 0.62 | 1.60
anisoAngleRadians | $\phi_{A}$ | 0.65 | 0.53 | 0.73 | 0.58 | 0.73
aniso1 | $\gamma_{2}$ | 0.70 | 0.30 | 1.28 | NA | NA
aniso2 | $\gamma_{3}$ | 2.57 | 1.55 | 3.56 | NA | NA
boxcox | $\lambda$ | 0.50 | 0.34 | 0.66 | 0.34 | 0.65

##  7 Example: Soil mercury data

To illustrate our method with a larger data set, we apply the method on the
considerably larger soil mercury dataset which has 829 observations and 29
predictors. The soil mercury data is obtained from “The EuroGeoSurveys -
FOREGS Geochemical Baseline Database" [15], available for download at
http://weppi.gtk.fi/publ/foregsatlas/ForegsData.php. The mercury concentration
measurement in topsoils is plotted in Figure 10(a) The covariate $X$ contains:
elevation, land category (see Figure 10(b) and 10(c)), night-time lights,
vegetation index.

(a) Mercury in soil (b) Land type (c) Land type

Figure 10: Soil mercury measurement with locations in Europe.

The procedure is identical to that used in the Swiss rainfall data example. We
fit 4 models using geostatsp::lgm() to obtain the MLEs, the first using the
MLE of the Matérn shape parameter and final three fixing the parameter at 0.5,
0.8, and 1. gpuLik::configParams() configured the candidate points based on
the MLEs of each model. The contour levels selected are the same as in the
Swiss rainfall data example, the parameter configuration matrix has 12316 rows
(including MLEs from the 4 fitted LGMs, $12\cdot 726$ configurations based on
the hgRes1 object below, and $3\cdot 10\cdot 120$ configurations based on the
other 3 models) and 5 columns. The values for boxcox is configured by choosing
31 equally-spaced values between the 0.01th (b[1]) and 0.99th (b[2]) quantiles
of the approximated Gaussian distribution. We obtained the log-likelihoods by
gpuLik::likfitLgmGpu() with the code below.

    R> hgRes1 = geostatsp::lgm(formula = HG ~ elevation + land + night + evi,
    +      data = hgm, grid = 20, covariates = covList, fixBoxcox = FALSE,
    +      fixShape = FALSE, fixNugget = FALSE, reml = FALSE, aniso = TRUE)

    R> reParamaters <- gpuLik::configParams(list(hgRes1, hgRes2, hgRes3, hgRes4),
    +      alpha = c(0.00001, 0.01, 0.1, 0.2, 0.25, 0.3, 0.5, 0.8, 0.9, 0.95,
    +          0.99, 0.999))
    R> paramsUse <- reParamaters$representativeParamaters[, 1:5]
    R> b <- reParamaters$boxcox

    R> result <- gpuLik::likfitLgmGpu(model = hgRes1, params = paramsUse,
    +      boxcox = seq(b[1], b[9], len = 31), cilevel = 0.95, type = "double",
    +      convexHullForBetas = FALSE, NparamPerIter = 400, Nglobal = c(256,
    +          256), Nlocal = c(16, 16), NlocalCache = 2800, verbose = c(1,
    +          0))

(a) $\sqrt{\phi_{X}\phi_{Y}}$, $\gamma_{1}$ (b) $\kappa$ (c) $\nu^{2}$ (d)
$\gamma_{2}$ (e) $\gamma_{3}$ (f) anisotropy ratio, $\phi_{R}$ (g) anisotropy
radians, $\phi_{A}$ (h) Legend

Figure 11: Profile log-likelihoods for the correlation parameters estimated
from the Soil mercury data. Vertical axes are log likelihood minus the
likelihood at the MLE. Dashed lines show the MLE and 95% confidence intervals.
Point clouds are representative points and their likelihoods, with colours
showing the $\alpha$ quantile of the contour each point was drawn from.

(a) $(\gamma_{2},\gamma_{3})$ (b) $(\phi_{R},\phi_{A})$ (c)
$(\gamma_{1},\nu^{2})$ (d) $(\gamma_{1},\gamma_{2})$ (e)
$(\gamma_{1},\gamma_{3})$ (f) $(\gamma_{1},\log(\kappa))$ (g)
$(\nu^{2},\gamma_{2})$ (h) $(\nu^{2},\gamma_{3})$

Figure 12: Two-dimensional profile log-likelihood contour plots for pairs of
model parameters, computed from the Soil mercury data. Circular plotting
symbols show locations of representative points on convex hull.

The soil mercury data is near isotropic, and the CI for the anisotropy radians
parameter covers the entire domain $(-\pi/2,\pi/2)$. The 2-dimensional profile
log-likelihood contour plots (a), (b) in Figure 12, using the internal
parameters does make the quadratic approximation for the likelihoods much
better. Parameter estimates and CIs are compared to the Wald CI’s in Table 5.
For most of the parameters, profile likelihood gives wider CIs than Wald does,
for parameters such as anisotropy radians, anisotropy ratio and shape,
_gpuLik_ gives better CIs as it takes into account the parameter boundaries.
Again, very similar CI’s for the transformation parameter are obtained by the
two methods.

(a) Standard deviation, $\sigma$ (b) Intercept, ($\beta_{0}$) (c)
Elevation,$\beta_{1}*10000$ (d) Box-Cox, $\lambda$ (e) night, $\beta_{16}$ (f)
landMosaic_grassland,$\beta_{5}$

Figure 13: Profile log-likelihoods for the standard deviation, some of the regression parameters, and Box Cox transformation parameter computed from the Soil mercury data. Horizontal dashed lines on each plot define approximate 95% confidence intervals and the MLE.  Table 5: 95 % CI’s for model parameters in the LGM fit to the Soil mercury data. Shown are the GPU-computed profile likelihood CI’s from _gpuLik_ and the traditional Wald-based CI’s computed by _geostatsp_. | likelihood-based | Wald
---|---|---
| notation | estimate | ci0.025 | ci0.975 | ci0.025 | ci0.975
(Intercept) | $\beta_{1}$ | -6.63 | -8.90 | -4.98 | -7.80 | -5.45
elevation*10000 | $\beta_{2}*10000$ | 2.72 | -1.02 | 6.57 | -1.00 | 6.40
broadleaved deciduous forest | $\beta_{3}$ | -0.28 | -0.65 | 0.08 | -0.63 | 0.08
Mosaic cropland / vegetation | $\beta_{4}$ | -0.08 | -0.43 | 0.25 | -0.42 | 0.25
Open needleleaved deciduous | $\beta_{5}$ | -0.63 | -1.21 | -0.10 | -1.16 | -0.10
Mosaic grassland | $\beta_{6}$ | -0.08 | -0.50 | 0.33 | -0.49 | 0.33
Mosaic vegetation / cropland | $\beta_{7}$ | -0.37 | -0.80 | 0.03 | -0.78 | 0.03
mixed broadleaved | $\beta_{8}$ | -0.51 | -1.01 | -0.04 | -0.97 | -0.05
Sparse vegetation | $\beta_{9}$ | -1.02 | -1.60 | -0.50 | -1.52 | -0.52
herbaceous vegetation | $\beta_{10}$ | 0.00 | -0.58 | 0.59 | -0.57 | 0.58
Mosaic forest or shrubland | $\beta_{11}$ | -0.28 | -0.91 | 0.32 | -0.88 | 0.31
needleleaved evergreen forest | $\beta_{12}$ | -0.01 | -0.60 | 0.57 | -0.58 | 0.56
shrubland | $\beta_{13}$ | 0.20 | -0.49 | 0.90 | -0.48 | 0.88
grassland or woody vegetation | $\beta_{14}$ | -0.69 | -1.64 | 0.20 | -1.58 | 0.20
Artificial surfaces | $\beta_{15}$ | -0.44 | -2.38 | 1.47 | -2.32 | 1.45
Water bodies | $\beta_{16}$ | -1.17 | -3.78 | 1.35 | -3.67 | 1.34
night*1000 | $\beta_{17}*1000$ | 1.57 | -0.34 | 3.57 | -0.32 | 3.44
evi | $\beta_{18}$ | 2.07 | 1.00 | 3.28 | 1.03 | 3.11
sdSpatial | $\sigma$ | 1.48 | 1.01 | 2.45 | 0.91 | 2.40
range/1000 | $\phi_{X}/1000$ | 1212.73 | 450.28 | 22100.81 | 220.56 | 6667.96
combinedRange/1000 | $\sqrt{\phi_{X}\phi_{Y}}/1000$ | 1138.60 | 365.27 | 19530.13 | NA | NA
anisoRatio | $\phi_{R}$ | 1.13 | 1.00 | 2.08 | 0.64 | 2.00
shape | $\kappa$ | 0.20 | 0.07 | 0.75 | -0.14 | 0.54
nugget | $\nu^{2}$ | 0.37 | 0.00 | 1.17 | NA | NA
sdNugget | $\tau$ | 0.90 | 0.00 | 1.60 | 0.38 | 2.13
anisoAngleRadians | $\phi_{A}$ | -0.82 | -1.57 | 1.57 | -2.55 | 0.92
aniso1 | $\gamma_{2}$ | -0.02 | -0.73 | 0.71 | NA | NA
aniso2 | $\gamma_{3}$ | -0.37 | -1.00 | 0.64 | NA | NA
boxcox | $\lambda$ | -0.23 | -0.29 | -0.18 | -0.29 | -0.17

##  8 Discussion

Inference for anisotropic geostatistical models is challenging due to the
large number of non-linear parameters involved (two range parameters,
anisotropic parameters, nugget effect), and more so when the Matérn shape
parameter and Box-Cox parameter are treated as unknown. Much of the software
for geostatistical analysis does not have the capability to compute profile
likelihoods for the model parameters, and instead uses plug-in information-
based (Wald) CIs for regression parameters (with the correlation parameters
treated fixed at their MLEs). These CIs do not take into account how
uncertainty in other parameters’ affects the uncertainty of $\beta$’s.
Furthermore, information-based CIs are not always possible, since the
information matrix is not positive definite when for example, the parameter
estimate is close to a boundary.

Computing and showing profile likelihoods can be helpful for evaluating
spatial model parameters, and likelihood-based CIs would provide wider and
better coverage than information-based CIs. Obtaining profile likelihood for
any of the parameters requires a large amount of repeated evaluation of the
likelihoods, this computation-heavy work is more easily tackled by parallel
computing with a GPU. GPUs are designed to process thousands of parallel
threads simultaneously, and are good at the repetitive likelihood evaluation
and can considerably reduce the run-time. The GPU-based technique introduced
here could be applied to other transformed-Gaussian models with many variance
parameters, such as hierarchical models with multiple random effects,
smoothing models with latent random walks, or more complex spatial models with
additional random effects.

Our work focused on using dense matrix spatial algorithms without
approximations, leveraging GPU’s to enable the standard Cholesky-based
likelihood evaluation to be performed on the reasonably large soil mercury
dataset with 800 spatial locations. The main computational limiting factor for
scaling up to larger datasets is GPU memory. Storing 400 variance matrices in
double precision for 10,000 spatial locations would require 300GB of memory,
which at the time of writing is only available on specialized and costly GPU
hardware. A second constraint would be the time taken to maximize the
likelihood function, which requires multiple sequential evaluations of the
likelihood on the CPU. Thirdly, large datasets with many points close to one
another (relative to the range parameter) could have variance matrices which,
after rounding errors, are not positive definite. However, any dataset for
which the MLE’s can be obtained using the existing dense matrix Cholesky
methodology could be accommodated with this GPU profile likelihood approach,
provided a sufficiently powerful GPU is used.

When datasets are larger (or time constraints more severe) than a dense
Cholesky algorithm can accommodate, approximations such as the SPDE approach
[16], spectral approximations [17], low rank approximations [18, 19], or
covariance tapering [20, 21] are necessary. The statistical methodology of
creating representative points and evaluating likelihoods in parallel is
perfectly compatible with such approximations, however implementing them would
require new GPU code designed to carry out the parallel computations
efficiently. The SPDE approach has become the dominant approximation for use
with Bayesian inference, and a sparse matrix equivalent to the batch matrix
functions in the current _gpuLik_ package would enable GPU profile likelihoods
for larger datasets using the SPDE model.

Extending the methodology to latent-Gaussian models, with a non-Gaussian
response distribution, would be a useful but more complex undertaking. There
is no closed-form expression for the likelihood and a Laplace approximation as
used by [22] and [10] would be the simplest approach. Achieving this on a GPU
would require an efficient batch implementation of ’inner’ optimization of the
joint probability $\pi(Y,U)$, which would be a substantial undertaking.

In summary, our method leveraging parallel likelihood evaluations on the GPU
is a feasible and convenient method for understanding the uncertainty in the
covariance parameters of anisotropic spatial models and propagating this
uncertainty through to inference on the $\beta$ regression parameters. This
includes the often disregarded Matérn shape parameter, we found that although
the profile likelihood of $\kappa$ is often flat, it generally has an
identifiable maximum and can be profiled out. Further, we have shown that
frequentist inference using exact form of the likelihood (without
approximating the variance matrix) is feasible for moderately large spatial
datasets.

##  9 Appendix: Restricted maximum profile likelihoods

One problem of the maximum likelihood estimation is that variance is
underestimated. Think of the simple case when the covariance matrix of $Y$ is
$\sigma^{2}I$, then

$\displaystyle\hat{\sigma}^{2}=\frac{1}{n}(y-X\hat{\beta})^{{{}^{{}_{\sf
T}}}}(y-X\hat{\beta}),$
$\displaystyle\mathbb{E}[\hat{\sigma}^{2}]=\frac{n-p}{n}\sigma^{2}<\sigma^{2},$

where $p=\text{rank}(X)$, is the number of elements in $\beta$. The estimation
bias comes from not taking into account the degree of freedom used for
estimating $\beta$. Restricted maximum likelihood (REML) developed by [23] and
[24] is an approach that produces unbiased estimators or less biased
estimators than MLE in general. REML eliminates the influence of $X$ on
$\hat{\sigma}^{2}$ and $\hat{\omega}$. In this method, we find a matrix $A$ of
full rank such that $AX=0$. The $n\times n$ projection matrix
$S=I-X(X^{{{}^{{}_{\sf T}}}}X)^{-1}X^{{{}^{{}_{\sf T}}}}$ satisfies $SX=0$,
however, $S$ is degenerate as $\text{rank}(S)=n-p$, so we take just $n-p$
linearly independent rows from $S$ to make the matrix $A_{(n-p)\times n}$.
Transform the observations $Y$ linearly to $Y^{*}=AY$, then
$\mathbb{E}(Y^{*})=AX\beta=0$. Recall in (1) we obtained $Y\sim
N_{n}(X\beta,\sigma^{2}V),$ thus

$Y^{*}\sim N_{n}(0,\sigma^{2}AVA^{{{}^{{}_{\sf T}}}}).$

The principle of the REML method is to estimate the variance parameter
$\sigma^{2}$ by maximizing the restricted log-likelihood function given below,
where $y^{*}$ represents a realization of the random variable $Y^{*}$,

$-2\ell^{*}(\sigma^{2},\omega;y^{*})=y^{*}{{{}^{{}_{\sf
T}}}}(\sigma^{2}AVA^{{{}^{{}_{\sf
T}}}})^{-1}y^{*}+\log|\sigma^{2}AVA^{{{}^{{}_{\sf T}}}}|+(n-p)\log(2\pi),$
(14)

Differentiating (14) with regard to $\sigma^{2}$ and setting it equal to zero
gives the unbiased REML estimate for $\sigma^{2}$

$\displaystyle\hat{\sigma}_{reml}^{2}(\omega)=\frac{y^{*}{{{}^{{}_{\sf
T}}}}(AVA^{{{}^{{}_{\sf T}}}})^{-1}y^{*}}{n-p}.$

The REML criterion is based on the likelihood function of $y^{*}$ (14), in
which $X$ does not appear. REML bypasses estimating $\beta$ and can therefore,
produce unbiased estimates for $\sigma^{2}$.

We can write (14) in terms of the observed data $y$ using the two identities
derived by [25] in Section 5.1:

$|A(\sigma^{2}V)A^{{{}^{{}_{\sf T}}}}|=|\sigma^{2}V||X^{{{}^{{}_{\sf
T}}}}(\sigma^{2}V)^{-1}X|$

and

$y^{*}(\sigma^{2}AVA^{{{}^{{}_{\sf
T}}}})^{-1}y^{*}=(y-X\hat{\beta})^{{{}^{{}_{\sf
T}}}}(\sigma^{2}V)^{-1}(y-X\hat{\beta}),$

where $\hat{\beta}=(X^{{{}^{{}_{\sf T}}}}V^{-1}X)^{-1}X^{{{}^{{}_{\sf
T}}}}V^{-1}y$. Then the restricted log-likelihood of $y$ is

$\displaystyle-2\ell^{*}(\sigma^{2},\omega;y)$
$\displaystyle=(y-X\hat{\beta})^{{{}^{{}_{\sf
T}}}}(\sigma^{2}V)^{-1}(y-X\hat{\beta})+(n-p)\log\sigma^{2}+\log|V|$
$\displaystyle+\log|X^{{{}^{{}_{\sf T}}}}V^{-1}X|+n\log(2\pi).$

If Box-Cox transformation is applied on $y$, then the general restricted log-
likelihood of $y$ incorporating the Box-Cox transformation parameter would be

$\displaystyle-2\ell^{*}(\sigma^{2},\omega,\lambda;y)$
$\displaystyle=(y^{\prime}-X\hat{\beta})^{{{}^{{}_{\sf
T}}}}(\sigma^{2}V)^{-1}(y^{\prime}-X\hat{\beta})+(n-p)\log\sigma^{2}+\log|V|$
$\displaystyle+\log|X^{{{}^{{}_{\sf
T}}}}V^{-1}X|+n\log(2\pi)-2(\lambda-1)\sum_{i=1}^{n}\log y_{i},$ (15)

which does not depend on the choice of matrix $A$. Differentiate (15) with
respect to $\sigma^{2}$ and set it equal to zero, we have the MLE
$\hat{\sigma}_{reml}^{2}$ in terms of the transformed data $y^{\prime}$ given
parameters $\omega$ and $\lambda$,

$\hat{\sigma}_{reml}^{2}(\omega,\lambda)=\frac{(y^{\prime}-X\hat{\beta})^{{{}^{{}_{\sf
T}}}}V^{-1}(y^{\prime}-X\hat{\beta})}{n-p}.$ (16)

Substitute (16) into (15), leading to the profile restricted log-likelihood
for $(\omega,\lambda)$

$\displaystyle-2\ell^{*}_{p}(\omega,\lambda;y)$
$\displaystyle=(n-p)\log\frac{(y^{\prime}-X\hat{\beta})^{{{}^{{}_{\sf
T}}}}V^{-1}(y^{\prime}-X\hat{\beta})}{n-p}+\log|V|+\log|X^{{{}^{{}_{\sf
T}}}}V^{-1}X|$ $\displaystyle-2(\lambda-1)\sum_{i=1}^{n}\log
y_{i}+n\log(2\pi)+n-p.$ (17)

Compared to the “ml” PLL given in (2.2), the effective dimension is reduced
from $n$ of $\ell(\cdot)$ to $(n-p)$ of $\ell^{*}_{p}(\cdot)$ , this
distinction is important if $p$ is large.

###  9.1 Profile likelihood for $\beta_{p}$ using REML

If REML method is used, then
$\ell_{p}(\beta_{i},\hat{\tilde{\beta}}_{\beta_{i}},\hat{\sigma}_{\beta_{i}},\omega,\lambda;y)$
would be equivalent to

$\displaystyle-2\ell^{*}_{p}(\omega,\lambda;y)$
$\displaystyle=(n-p)\log\frac{(\tilde{y}-\tilde{X}\hat{\tilde{\beta}}_{\omega,\lambda})^{{{}^{{}_{\sf
T}}}}V^{-1}(\tilde{y}-\tilde{X}\hat{\tilde{\beta}}_{\omega,\lambda})}{n-p}+\log|V|+\log|X^{{{}^{{}_{\sf
T}}}}V^{-1}X|$ $\displaystyle-2(\lambda-1)\sum_{i=1}^{n}\log
y_{i}+n\log(2\pi)+n-p.$ (18)

## References

  * [1] George EP Box and David R Cox.  An analysis of transformations.  Journal of the Royal Statistical Society: Series B (Methodological), 26(2):211–243, 1964.
  * [2] Dionissios T Hristopulos.  Random fields for spatial data modeling: a primer for scientists and engineers.  Springer Nature, 2020.
  * [3] Michael L Stein.  Interpolation of spatial data: some theory for kriging.  Springer Science & Business Media, 1999.
  * [4] Peter J Diggle and Emanuele Giorgi.  Model-based geostatistics for global public health: methods and applications.  Chapman and Hall/CRC, 2019.
  * [5] Kathryn A Haskard, Brian R Cullis, and Arūnas P Verbyla.  Anisotropic matérn correlation and spatial prediction using reml.  Journal of agricultural, biological, and environmental statistics, 12(2):147–160, 2007.
  * [6] Charlotte Baey, Paul-Henry Cournède, and Estelle Kuhn.  Asymptotic distribution of likelihood ratio test statistics for variance components in nonlinear mixed effects models.  Computational Statistics & Data Analysis, 135:107–122, 2019.
  * [7] Peter M Visscher.  A note on the asymptotic distribution of likelihood ratio tests to test variance components.  Twin research and human genetics, 9(4):490–495, 2006.
  * [8] Patrick E Brown.  Model-based geostatistics the easy way.  Journal of Statistical Software, 63:1–24, 2015.
  * [9] PD Dauncey, M Kenzie, N Wardle, and GJ Davies.  Handling uncertainties in background shapes: the discrete profiling method.  Journal of Instrumentation, 10(04):P04015, 2015.
  * [10] Håvard Rue, Sara Martino, and Nicolas Chopin.  Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations.  Journal of the royal statistical society: Series b (statistical methodology), 71(2):319–392, 2009.
  * [11] Patrick E Brown Kamal Rai.  Estimation and prediction in spatial models with block composite likelihoods.  2022\.  Manuscript in preparation.
  * [12] Hao Zhang.  Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics.  Journal of the American Statistical Association, 99(465):250–261, 2004.
  * [13] Richard Arnold Johnson, Dean W Wichern, et al.  Applied multivariate statistical analysis, volume 6.  Pearson London, UK:, 2014.
  * [14] Patrick E Brown.  Package ‘geostatsp’.  2021\.
  * [15] R Salminen, MJ Batista, M Bidovec, A Demetriades, B De Vivo, W De Vos, M Duris, A Gilucis, V Gregorauskiene, J Halamic, et al.  Geochemical atlas of europe. part 1—background information, methodology and maps. geological survey of finland.  Beschikbaar op: http://www. gtk. fi/publ/foregsatlas/index. php [30/03/2010], 2005.
  * [16] Finn Lindgren, Håvard Rue, and Johan Lindström.  An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(4):423–498, 2011.
  * [17] Benjamin M Taylor, Tilman M Davies, Barry S Rowlingson, and Peter J Diggle.  lgcp an r package for inference with spatio-temporal log-gaussian cox processes.  arXiv preprint arXiv:1110.6054, 2011.
  * [18] Sudipto Banerjee, Alan E Gelfand, Andrew O Finley, and Huiyan Sang.  Gaussian predictive process models for large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(4):825–848, 2008.
  * [19] Noel Cressie and Gardar Johannesson.  Fixed rank kriging for very large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(1):209–226, 2008.
  * [20] Reinhard Furrer, Marc G Genton, and Douglas Nychka.  Covariance tapering for interpolation of large spatial datasets.  Journal of Computational and Graphical Statistics, 15(3):502–523, 2006.
  * [21] Cari G Kaufman, Mark J Schervish, and Douglas W Nychka.  Covariance tapering for likelihood-based estimation in large spatial data sets.  Journal of the American Statistical Association, 103(484):1545–1555, 2008.
  * [22] Alex Stringer, Patrick E Brown, and Jamie Stafford.  Fast, scalable approximations to posterior distributions in extended latent gaussian models.  arXiv preprint arXiv:2103.07425, 2021.
  * [23] H Desmond Patterson and Robin Thompson.  Recovery of inter-block information when block sizes are unequal.  Biometrika, 58(3):545–554, 1971.
  * [24] David Roxbee Cox and Nancy Reid.  Parameter orthogonality and approximate conditional inference.  Journal of the Royal Statistical Society: Series B (Methodological), 49(1):1–18, 1987.
  * [25] SR Searle, RL Quaas, et al.  A notebook on variance components: A detailed description of recent methods of estimating variance components, with applications in animal breeding.  1978\.