7 Continuous Spatial Variation 7.1 Geostatistical models 7.3 Model Fitting Via The Variogram

7.2 Model Fitting Via Maximum Likelihood

The model in Equation 7.1 can be written in matrix form as:

Y=\mathrm{MVN}(Z\beta,\Sigma_{\sigma,\phi}+\tau^{2}I),

where $Z$ is the design matrix of covariate measurements, $\beta$ are the covariate effects, $\Sigma_{\sigma,\phi}$ is the variance of the spatial process $S$ at the locations $x_{1},\ldots,x_{n}$ , $I$ is the ( $n\times n$ ) identity matrix and $\tau^{2}$ is the variance of the spatially uncorrelated random effects. We can write $\Sigma_{\sigma,\phi}=\sigma^{2}R_{\phi}$ , where the $i, j$ -element of the matrix $R$ is formed by evaluating $\rho(||x_{i}-x_{j}||)$ and $||\,\cdot\,||$ denotes Euclidean distance. Put $\nu^{2}=\tau^{2}/\sigma^{2}$ and let $M_{\nu,\phi}=R_{\phi}+\nu^{2}I$ , we re-write the model above as

Y\sim\mathrm{MVN}(Z\beta,\sigma^{2}M_{\nu,\phi}).

(7.3)

The log-likelihood for a set of data $Y$ is then:

\log\pi(\text{data}|\beta,\sigma,\nu,\phi)=\text{const}-\frac{1}{2}\log|\sigma% ^{2}M_{\nu,\phi}|-\frac{1}{2\sigma^{2}}(Y-Z\beta)^{T}M_{\nu,\phi}^{-1}(Y-Z% \beta).

(7.4)

For given $\nu$ and $\phi$ , we can maximise the expression above exactly using estimated generalised least squares (EGLS); this yields estimates of $\beta$ and $\sigma$ , denoted respectively $\hat{\beta}^{(1)}$ and $\hat{\sigma}^{(1)}$ . Now, given $\hat{\beta}^{(1)}$ and $\hat{\sigma}^{(1)}$ , we can numerically maximise the profile log-likelihood to get estimates for $\nu$ and $\phi$ , denoted $\hat{\nu}^{(1)}$ and $\hat{\phi}^{(1)}$ . Next, given $\hat{\nu}^{(1)}$ and $\hat{\phi}^{(1)}$ , we can maximise exactly to obtain $\hat{\beta}^{(2)}$ and $\hat{\sigma}^{(2)}$ and again numerically maximise the resulting profile likelihood to obtain $\hat{\nu}^{(2)}$ and $\hat{\phi}^{(2)}$ .

Iterating this procedure, we find the MLEs of $\beta$ , $\sigma$ , $\nu$ and $\phi$ and the estimated Hessian, from which we can derive confidence intervals. There are two issues worth noting.

•

As with all complex maximum likelihood routines, ideally one should attempt the maximisation process from several starting values. Note that the variogram (see below) can be used to generate initial values.
•

Since the above algorithm requires the invertion of an $n\times n$ matrix, it can be prohibitively slow for large datasets

Regardless of which method of estimation is chosen, it is always a good idea to plot the variogram (see below) to check if the MLEs of $\beta$ , $\sigma$ , $\nu$ and $\phi$ look sensible.

Example 7.4.

Estimated Generalised Least Squares (I) Here we show how to derive the MLE for $\beta$ in a model with correlated residuals. Start with $Y=X\beta+\epsilon$ , where $\epsilon\sim\mathrm{N}(0,\Sigma)$ and $\Sigma$ is a covariance matrix. Since $\Sigma$ is symmetric and positive definite, we can use the eigendecomposition to write it as $\Sigma=U\Lambda U^{T}$ , which means $\Sigma^{-1}=U\Lambda^{-1}U^{T}$ . We write $U\Lambda^{-1}U^{T}=U\Lambda^{-1/2}(\Lambda^{-1/2})^{T}U^{T}=VV^{T}$ . Pre-multiplying the original equation by $V^{T}$ gives $V^{T}Y=V^{T}X\beta+V^{T}\epsilon$ , so $\mathbb{V}(V^{T}\epsilon)=V^{T}\Sigma V=I$ , where $I$ is the identity matrix. Therefore, in pre-multiplying by $V^{T}$ , we reduce our original model to one that can be solved by ordinary least squares and we can derive the EGLS estimate of $\beta$ by solving the matrix equation $V^{T}Y=V^{T}X\beta$ for $\beta$ . Premultiplying by $X^{T}V$ gives

$\displaystyle X^{T}VV^{T}Y$	$\displaystyle=$	$\displaystyle X^{T}VV^{T}X\beta$
$\displaystyle X^{T}\Sigma^{-1}Y$	$\displaystyle=$	$\displaystyle X^{T}\Sigma^{-1}X\beta$
$\displaystyle X^{T}\Sigma^{-1}Y$	$\displaystyle=$	$\displaystyle X^{T}\Sigma^{-1}X\beta$
$\displaystyle\hat{\beta}_{\text{EGLS}}$	$\displaystyle=$	$\displaystyle(X^{T}\Sigma^{-1}X)^{-1}X^{T}\Sigma^{-1}Y$	(7.5)

Example 7.5.

Estimated Generalised Least Squares (II) Here we show how to derive the MLE for $\beta$ and $\sigma^{2}$ in Equation 7.3. Applying EGLS (i.e. Equation 7.5) to this model yields $\hat{\beta}=(X^{T}(\sigma^{2}M_{\nu,\phi})^{-1}X)^{-1}X^{T}(\sigma^{2}M_{\nu,% \phi})^{-1}Y$ ; the $\sigma$ s cancel, hence

\hat{\beta}=(X^{T}M_{\nu,\phi}^{-1}X)^{-1}X^{T}M_{\nu,\phi}^{-1}Y.

Now suppose that we have obtained the MLE for $\beta$ so that conditional on this parameter, $Y\sim\mathrm{MVN}(Z\hat{\beta},\sigma^{2}M_{\nu,\phi})$

First observe that $Y-\mathbb{E}(Y)=Y-Z\hat{\beta}\sim\mathrm{MVN}(0,\sigma^{2}M_{\nu,\phi})$ and so $M_{\nu,\phi}^{-1/2}(Y-Z\hat{\beta})\sim\mathrm{MVN}(0,\sigma^{2})$ i.e. the variance of the column vector $v=M_{\nu,\phi}^{-1/2}(Y-Z\hat{\beta})$ is an estimate of $\sigma^{2}$ . Since $\mathbb{E}(v)=0$ we compute this as $\frac{1}{n}v^{T}v$ , that is,

\hat{\sigma}^{2}=\frac{1}{n}(Y-Z\hat{\beta})^{T}M_{\nu,\phi}^{-1}(Y-Z\hat{% \beta}).