\documentclass[a4paper]{article}

\usepackage{Sweave}

\usepackage[utf8]{inputenc}

\usepackage[figuresright]{rotating}
\usepackage{tabu}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{mathtools}
\usepackage{pifont}
\usepackage{lineno}
\usepackage{graphicx}
\usepackage{enumerate}
\usepackage{rotating}
\usepackage[round]{natbib}
\usepackage{scrextend}
\usepackage{datetime}
\usepackage{array}
\usepackage{multicol}
\usepackage[justification = raggedright]{subcaption}
%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\proglang}[1]{{\textsf{#1}}}
\newcommand{\R}{\proglang{R }}
\newcommand{\pkg}[1]{{\normalfont\fontseries{b}\selectfont #1}}
\newcommand{\pbs}[1]{\let\tmp\\#1\let\\\tmp} 
%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\s}{\mathbf{s}}
\newcommand{\h}{\mathbf{h}}
\newcommand{\ttt}{\mathbf{t}}
\newcommand{\0}{\mathbf{0}}
\newcommand{\A}{\mathbf{A}}
\newcommand{\G}{\mathbf{G}}
\newcommand{\COV}{\mathbf{COV}}
\newcommand{\VAR}{\mathbf{VAR}}
\newcommand{\D}{\mathcal{D}}
\newcommand{\E}{\mathbf{E}}
%\newcommand{\Re}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\x}{\mathbf{x}}
\newcommand{\bg}{\boldsymbol{\gamma}}
\newcommand{\bo}{\boldsymbol{\omega}}
\newcommand{\be}{\boldsymbol{\epsilon}}
\newcommand{\bL}{\mathbf{\Lambda}}
\newtheorem{definition}[subsubsection]{Definition}
\newtheorem{remark}[section]{Remark}
\numberwithin{equation}{section}
\newtheorem{definitionA}[subsubsection]{Definition}
%%%%%%%%%%%%%%%%%%%%%%%%%
\SweaveOpts{prefix.string=figures/fig}
\SweaveOpts{keep.source=FALSE} 
\SweaveOpts{highlight=TRUE, tidy=TRUE, keep.space=TRUE, keep.blank.space=FALSE, keep.comment=TRUE}
\SweaveOpts{eps=TRUE, pdf=FALSE, png=FALSE}

\defcitealias{lu2005test}{LZ}
\defcitealias{guan2004isotropy}{GSC}
\defcitealias{maity2012test}{MS}

% \VignetteIndexEntry{spTest-vig}
% \VignetteDepends{geoR, fields, splines}
% \VignetteKeyword{anisotropy}

\title{spTest: An R Package Implementing Nonparametric Tests of Isotropy}
\author{Zachary D. Weller}

\begin{document}
%\SweaveOpts{concordance=TRUE}


\maketitle

\begin{abstract}
An important step in modeling spatially-referenced data is appropriately specifying the second order properties of the random field. A scientist developing a model for spatial data has a number of options regarding the nature of the dependence between observations. One of these options is deciding whether or not the dependence between observations depends on direction, or, in other words, whether or not the spatial covariance function is isotropic. Isotropy implies that spatial dependence is a function of only the distance and not the direction of the spatial separation between sampling locations. A researcher may use graphical techniques, such as directional sample semivariograms, to determine whether an assumption of isotropy holds. These graphical diagnostics can be difficult to assess, subject to personal interpretation, and potentially misleading as they typically do not include a measure of uncertainty. In order to escape these issues, a hypothesis test of the assumption of isotropy may be more desirable. To avoid specification of the covariance function, a number of nonparametric tests of isotropy have been developed using both the spatial and spectral representations of random fields. Several of these nonparametric tests are implemented in the \R package \pkg{spTest}, available on \texttt{CRAN}. We demonstrate how graphical techniques and the hypothesis tests programmed in \pkg{spTest} can be used in practice to assess isotropy properties. A version of this vignette is also available in the Journal of Statistical Software \citep{spTestJSS}.
\end{abstract}

\section{Introduction}
An important step in modeling a spatial process is choosing the form of the covariance function. The form of the covariance function will have an effect on kriging as well as parameter estimates and the associated uncertainty \citep[pg.~127-135]{cressie1993}. A common simplifying assumption about the spatial covariance function is that it is isotropic, that is, the dependence between sampling locations depends only on the distance between locations and not on their relative orientation. This assumption may not always be reasonable; for example, wind may lead to directional dependence in environmental monitoring data. Misspecification of the second order properties can lead to misleading inference. \citet[pg.~87-90]{sherman2011spatial} and \citet{guan2004isotropy} summarize the effects of incorrectly specifying isotropy properties on kriging estimates through numerical examples. In order to choose an appropriate model, a statistician must first assess the nature of the spatial variation of his or her data. To check for anisotropy (directional dependence) in spatially-referenced data, a number of graphical techniques are available, such as directional sample variograms or rose diagrams \citep[pg.~149-154]{matheron1961precision,isaaks1989applied}. Spatial statisticians may have intuition about the interpretation and reliability of these diagnostics, but a less experienced user may desire evaluation of assumptions via a hypothesis test. 

A number of nonparametric tests of isotropy, which avoid the choice of a parametric covariance function, have been developed \citep[see, e.g.,][]{lu2001testing,guan2004isotropy, lu2005test, maity2012test}. In \citet{weller2015}, we provide a review of the different nonparametric methods available for testing isotropy and symmetry properties, including an extensive simulation study. In the current work, we aim to showcase the functionality of our \R \citep{rsoftware} package \pkg{spTest} \citep{spTest}, which implements several of the aforementioned methods. We use two real data examples to illustrate how the nonparametric hypothesis tests available in \pkg{spTest} can be used to assess isotropy properties in spatially-referenced data. The examples also demonstrate how graphical techniques and hypothesis tests can be used in a complementary role. The remainder of the paper is organized as follows: Section \ref{notation} establishes notation and definitions; Section \ref{sptest-functions} describes the functionality of the \pkg{spTest} package; Section \ref{applications} demonstrates how to use the functions in \pkg{spTest} in conjunction with graphical techniques on two different data sets; Section \ref{discussion} concludes the paper with a discussion.

\section{Notation}\label{notation}

\subsection{Spatial statistics}\label{spatial}
In Section~\ref{spatial} we review key definitions required for describing and performing nonparametric tests of isotropy. In particular we review the concepts of stationarity and isotropy in random spatial processes. Readers who are familiar with these concepts may wish to proceed to Section~\ref{nptests}. For additional background, see \citet{weller2015} or \citet{schab}. 

Let $\{Y(\s): \s \in \D \subseteq \Re^d\}$ be a random field (RF). For the remainder of the paper we assume the common case of $d = 2$. Let $\{\s_1, \ldots, \s_n\} \in \D$ be the finite set of locations at which the RF is observed, providing the random vector $\mathbf{Y} = (Y(\s_1), \ldots , Y(\s_n))^\top$. The sampling locations may follow one of several spatial sampling designs, for example, gridded locations, randomly and uniformly distributed locations, a cluster design, or any other general design. Note that there is a distinction between a lattice process and a geostatistical process observed on a grid \citep[e.g.,][pg.~6-10]{fuentes2010spectral, schab}, but we do not explore this distinction and will use the term grid throughout.

It is often of interest to infer the effect of covariates on the RF, deduce dependence structure, and/or predict $Y$ with associated uncertainty at new locations. To achieve these goals, we must specify the distributional properties of the spatial process. A common assumption is that the finite-dimensional joint distribution of $\{Y(\s): \s \in \D \subseteq \Re^2\}$ is multivariate normal, in which case we call the RF a Gaussian random field (GRF).

A RF is strictly stationary if its distribution is invariant under coordinate translation,
\begin{equation}\label{strict}
\mathbb{P}(Y(\s_1) < y_1, \cdots, Y(\s_k) < y_k ) = \mathbb{P}(Y(\s_1 + \h) < y_1, \cdots, Y(\s_k + \h) < y_k )
\end{equation}
for all spatial lags $\h = (h_1, h_2)^\top \in \Re^2$ and $k \in \mathbb{Z}^{+}$. A RF is weakly, or second order, stationary if
\begin{equation}\label{secondorder}
	\begin{aligned}
	\E(Y(\s)) &= \E(Y(\s + \h)) = \mu \text{ and} \\
	 \COV(Y(\s), Y(\s + \h)) &=  \COV(Y(\0), Y(\h)) = C(\h),
	\end{aligned}
\end{equation}	
where $C(\h) < \infty$ is called the covariance function or covariogram. Henceforth we assume $\mu = 0$. A RF that is weakly stationary \eqref{secondorder} has a constant mean, and the covariance between values at two spatial locations depends only on the spatial lag. Absolute coordinates are irrelevant for the mean, variance, and covariance of a weakly stationary RF. A RF is intrinsically stationary if 
\begin{equation} \label{intrinsic}
	\begin{aligned}
		\E(Y(\s + \h) - Y(\h)) &= 0 \text{ and }\\ 
		\VAR(Y(\s + \h) - Y(\s))  &= 2\bg(\h),
	\end{aligned}
\end{equation}	
where $2\bg(\h)$ is called the variogram function or variogram, and $\bg(\h) = \frac{1}{2}\VAR(Y(\s + \h) - Y(\s))$ is the semivariogram function or semivariogram. Note that intrinsic stationarity is defined in terms of the increments $Y(\s + \h) - Y(\s)$. 

The relationship between the different types of stationarity is given by
\begin{equation}\label{stationary_rel}
 \text{strict} \implies \text{weak} \implies \text{intrinsic}.
\end{equation}
In general the converse in \eqref{stationary_rel} is not true; however, the converse in \eqref{stationary_rel} holds under the assumption that the RF is a GRF. From \eqref{secondorder} and \eqref{intrinsic} it follows that $\bg(\h) = C(\0) - C(\h)$ when the RF is weakly stationary. This relationship implies that the second order properties of a weakly stationary RF can be viewed from the perspective of either the (semi)variogram or covariance function.

A common simplifying assumption about the spatial dependence structure is that it is isotropic.
\begin{definition}\label{iso}
A weakly stationary RF is isotropic if the covariance function, $C(\h)$, (or equivalently, the semivariogram function) of the RF depends on the lag vector $\h$ only through its Euclidean length, $||\h||$, i.e., $C(\h) = C_{0}(||\h||)$ for some function $C_{0}(\cdot)$ of a univariate argument.
\end{definition}
Isotropy implies that the dependence between any two observations depends only on the Euclidean distance between their sampling locations and not on their relative orientation. A spatial process that is not isotropic is called anisotropic. Figure~\ref{simulated_data} displays data simulated from a two stationary RFs, one which is isotropic and one that is anisotropic. The methods in the \R package \pkg{spTest} are designed to assist in determining whether or not an assumption of isotropy is reasonable.

<<label = settings, echo = FALSE, include = FALSE>>=
options(prompt = "R> ", continue = "+  ", width = 70)
@

<<label = libs1, echo = TRUE>>=
library("mvtnorm")
library("fields")
@

<<label=isotropicsim,fig=TRUE, include=FALSE>>=

matern.2d = function(h, phi, alpha, nu)
{

	a <- (pi*phi)/(2^(nu - 1)*gamma(nu + 1) * alpha^(2*nu))
	d <- h*alpha
	d[d == 0] <- 1e-10
	b <- (d)^nu
	c <- besselK(d, nu)
	
	cov.val <- a*b*c
	return(cov.val)
}

#Function to transform coordinates to create anisotropy
coords.aniso = function (coords, aniso.pars, reverse = FALSE) 
{
    coords <- as.matrix(coords)
    n <- nrow(coords)
    if (length(aniso.pars) != 2) 
        stop("argument aniso.pars must be a vector with 2 elements: the anisotropy angle and anisotropy ratio, respectively")
    psiA <- aniso.pars[1]
    psiR <- aniso.pars[2]
    if (psiR < 1) {
        psiR <- round(psiR, digits = 8)
        if (psiR < 1) 
            stop("anisotropy ratio must be greater than 1")
    }
    rm <- matrix(c(cos(psiA), -sin(psiA), sin(psiA), cos(psiA)), 
        ncol = 2)
    tm <- diag(c(1, 1/psiR))
    if (reverse) 
        coords.mod <- coords %*% solve(rm %*% tm)
    else coords.mod <- coords %*% rm %*% tm
    return(coords.mod)
}


set.seed(2017)

x = seq(0, 12, by = 0.25)
y = seq(0, 12, by = 0.25)

coords = expand.grid(x,y)
n = dim(coords)[1]

aniso.angle <- 3*pi/4
aniso.ratio <- 2
coordsA <- coords.aniso(coords, c(aniso.angle, aniso.ratio))
Da <-  as.matrix(dist(coordsA))

sigma.sq <- 1.0
alpha <- 1/2
sm <- 1
tau.sq <- 0

d <- seq(0, 20, by = 0.1)
phi <- sigma.sq * alpha^(2*sm) * gamma(sm + 1)/ (pi * gamma(sm))
mcov2 <- matern.2d(d, phi, alpha, sm)
er.index <- which( round(mcov2, 2) == 0.05)
er <- d[er.index]
er <- er[1]

D <- as.matrix(dist(coords))
R <- matern.2d(D, phi, alpha, sm)
R <- R + diag(tau.sq, nrow = n, ncol = n)

Ra <- matern.2d(Da, phi, alpha, sm)
Ra <- Ra + diag(tau.sq, nrow = n, ncol = n)

z <- rmvnorm(1,rep(0,n), R, method = c("chol"))
z <- z - mean(z)
zmat <- matrix(z, nrow = length(y), ncol = length(x), byrow = T)
zmat.iso <- zmat

x.shift = -102
y.shift = 37
par(mai = c(0.8, 0.4, 0.8, 0.2), pty = "s")
image(x+x.shift, y+y.shift, zmat,
ylab = "", xlab = "", col = two.colors(n = 256, start = "blue3", end = "red3", middle = "gray60"),
cex.main = 1.80, pty = "s")
map("state", add = T, lwd = 2)
mtext("Stationary and Isotropic Spatial Process", side = 3, line = 0.5, cex = 1.25, font = 2)
mtext("Longitude", side = 1, line = 2)
mtext("Latitude", side = 2, line = 2)
@

<<label=anisotropicsim,fig=TRUE, include=FALSE>>=

set.seed(2)
z <- rmvnorm(1,rep(0,n), Ra, method = c("chol"))
z <- z - mean(z)
zmat <- matrix(z, nrow = length(y), ncol = length(x), byrow = T)

zmat.anis <- zmat
zmat.iso <- zmat
x.shift <- -102
y.shift <-  37
par(mai = c(0.8, 0.4, 0.8, 0.2), pty = "s")
image(x+x.shift, y+y.shift, zmat,
ylab = "", xlab = "", col = two.colors(n = 256, start = "blue3", end = "red3", middle = "gray60"),
cex.main = 1.80, pty = "s")
map("state", add = T, lwd = 2)
mtext("Stationary and Anisotropic Spatial Process", side = 3, line = 0.5, cex = 1.25, font = 2)
mtext("Longitude", side = 1, line = 2)
mtext("Latitude", side = 2, line = 2)
@

\begin{figure}
	\begin{subfigure}[t]{0.5\linewidth}
	\centering
	\includegraphics[width = 3 in]{figures/fig-isotropicsim.eps}
	\caption{}
	\label{isotropic-fig}
	\end{subfigure}
	\begin{subfigure}[t]{0.5\linewidth}
	\centering
	\includegraphics[width = 3 in]{figures/fig-anisotropicsim.eps}
	\caption{}
	\label{anisotropic-fig}
	\end{subfigure}
\caption{Heat maps displaying simulated data from two stationary RFs. Figure~\ref{isotropic-fig} displays data simulated from an isotropic covariance function. Figure~\ref{anisotropic-fig} shows data simulated from an anisotropic covariance function. The existence of anisotropy in Figure~\ref{anisotropic-fig} is evidenced by bands of similar values running the northwest to southeast directions. These bands indicate that spatial dependence is stronger in the NW-SE direction than the SW-NE direction.}
\label{simulated_data}
\end{figure}

When modeling a RF, a typical assumption is that the spatial covariance, $C(\h))$, or semivariogram function, $\bg(\h)$, can be described by a parametric model. A number of methods are available for estimating the parameters of these models, for example, maximum likelihood or least squares \citep[pg.~90-97]{cressie1993}. Parametric models ensure that the covariance function is valid and provide parameters that can be interpreted as describing characteristics of the random field, such as the range of dependence or the direction of anisotropy \citep[pg.~141-152]{schab}. The methods in \pkg{spTest} are nonparametric tests of isotropy, which circumvent the choice of a parametric form for the covariance (semivariogram) function. A nonparametric test of isotropy escapes potential misspecification of a parametric covariance function and avoids the potential for having to estimate the covariance function twice (e.g., under the null and alternative hypothesis for a likelihood-ratio test).

A nonparametric test of isotropy requires a nonparametric estimator of second order properties of the RF. We discuss nonparametric estimation of the semivariogram function here and note that similar techniques can be used for nonparametric estimation of the covariance function. The classical moment-based estimator of the semivariogram \citep{matheron1962traite} is the sample semivariogram given by
	\begin{equation}\label{classical}
		\widehat{\gamma}(\h) = \frac{1}{2|\D(\h)|} \sum [Y(\s) - Y(\s + \h)]^2,
	\end{equation}
where the sum is over $\D(\h) = \{\s:\s, \s + \h \in \D\}$ and $|\D(\h)|$ is the number of elements in $\D(\h)$. The set $\D(\h)$ is the set of sampling location pairs that are separated by spatial lag $\h$. There are two important modifications to the estimator in \eqref{classical} that are pertinent to the methods described in this paper. First, for non-gridded sampling locations, the estimator needs to be modified to account for the fact that very few or no pairs of locations will be separated by a specific spatial lag, $\h$. One solution to this challenge is to specify a distance tolerance, $\epsilon$, such that lags having length $||\h|| \pm \epsilon$ are included in estimating the semivariogram at lag $\h$. Second, directional sample semivariograms can be estimated by using only observations that are separated by spatial lags in a specific direction. For example, to investigate potential anisotropy, we can compare sample semivariograms between the horizontal and vertical directions. For non-gridded sampling locations, very few pairs of locations will lie at a specific distance and directional lag, so we need to allow for both a distance and a directional tolerance when estimating the semivariogram. A common method for doing this is by using a product kernel smoother that smoothes over both the horizontal ($h_1$) and vertical ($h_2$) components of the spatial lag $\h = (h_1, h_2)^\top$.

Spatial RFs and their second order properties can also be expressed in the spectral (or frequency) domain using Fourier transforms. The spectral representation of RFs and their second order properties provides alternative methods for testing second order properties. For brevity we focus only on the methods in the spatial domain and refer the interested reader to \citet{weller2015} or \citet{fuentes2010spectral}. We note that that, in addition to the methods from the spatial domain, the nonparametric spectral methods from \citet{lu2005test} are also implemented in \pkg{spTest}.

\subsection{Nonparametric tests of isotropy}\label{nptests}
\citet{lu1994distributions} and \citet{lu2001testing} pioneered a popular approach to testing second-order properties when they used the joint asymptotic normality of the sample semivariogram computed at different spatial lags to evaluate symmetry and isotropy properties. The subsequent works of \citet{guan2004isotropy, guan2007asymptotic} and \citet{maity2012test} built upon these ideas and are the primary methods programmed in \pkg{spTest}. Here we give an overview of the tests in \citet{guan2004isotropy} and \citet{maity2012test}.

Under the null hypothesis that the RF is isotropic, it follows that the values of $\bg(\cdot)$ evaluated at any two spatial lags that have the same norm are equal, independent of the direction of the lags. For example, under the assumption of isotropy, $\bg((-6,0)) = \bg\left(\left(\sqrt{3}, \sqrt{3} \hspace{1 mm} \right)\right)$. To completely specify the null hypothesis of isotropy, theoretically, one would need to compare semivariogram values for an infinite set of lags. In practice, a small number of lags are specified. Then it is possible to test a hypothesis consisting of a set of linear contrasts of the form 
\begin{equation}\label{linearhypothesis}
H_0: \A\bg(\cdot) = \mathbf{0}
\end{equation}
as a proxy for the null hypothesis of isotropy, where $\A$ is a full row rank matrix \citep{lu2001testing}. For example, a set of lags, denoted $\bL$, commonly used in practice on gridded sampling locations with unit spacing is
\begin{equation}\label{lagset}
 \bL = \{ \h_1 = (1,0), \h_2 = (0,1), \h_3 = (1,1),  \h_4 = (-1,1) \}, 
 \end{equation}
 and the corresponding $\A$ matrix under $H_0: \A\bg(\bL) = \mathbf{0}$ is
\begin{equation}\label{amat}
	\A = \begin{bmatrix}
		1 & -1 & 0 & 0 \\
		0 & 0 & 1 & -1 
		\end{bmatrix}.
\end{equation}

In other words, using \eqref{lagset} and \eqref{amat} under the hypothesis \eqref{linearhypothesis}, we contrast the semivariogram values at lags $\h_1 = (1,0)$ and $\h_2 = (0,1)$, and likewise, contrast semivariogram values at lags $\h_3 = (1,1)$ and $\h_4 = (-1,1)$. An important step in detecting anisotropy is the choice of lags, $\bL$, as the test results will only hold for the set of lags considered \citep{guan2004isotropy}. While this choice is somewhat subjective, there are several considerations for determining the set of lags. First, in terms of Euclidean distance, short lags should be chosen. Short lags are preferred because estimates of the semivariogram at long lags tend to be more variable. Second, pairs of orthogonal lags should be contrasted because, for an anisotropic process, the directions of strongest and weakest spatial correlation will typically be orthogonal. For other considerations and more detailed guidelines regarding the choice of lags, see \citet{weller2015}, \citet{lu2001testing}, and \citet{guan2004isotropy}.

\begin{table}[h]
\begin{center}
\begin{tabular}{| l | | c | c | >{\centering}p{2.45 cm} | c | } 
\hline 
& \multicolumn{4}{c|}{\textbf{Hypothesis test properties}}   \\
\hline
 \textbf{Method} & \pkg{spTest} func'n & Design & Estimator & Primary user choices    \\ 
 \hline
  \citetalias{guan2004isotropy} & \texttt{GuanTestGrid} & grid & sample semivariogram  & lag set, window size    \\ 
 \hline
 \citetalias{lu2005test} & \texttt{LuTest} & grid   & periodogram &  $-$  \\ 
 \hline
 \citetalias{guan2004isotropy} & \texttt{GuanTestUnif} &  uniform  & kernel sample semivariogram  & lag set, window size, bandwidth   \\ 
 \hline
 \citetalias{maity2012test} &  \texttt{MaityTest} & any & kernel sample covariogram & lag set, block size   \\
\hline
\end{tabular}

\caption{Nonparametric tests of isotropy available in the \R package \pkg{spTest}, including \citetalias{guan2004isotropy}: \citet{guan2004isotropy}, \citetalias{lu2005test}: \citet{lu2005test}, and  \citetalias{maity2012test}: \citet{maity2012test}. The column ``\pkg{spTest} func'n" lists the name of the function used to implement the test in the \pkg{spTest} package. ``Design" indicates the spatial sampling design for which the test is valid. ``Estimator" describes the method used to estimate second order properties. Note that a ``uniform" design implies sampling locations must be uniformly distributed on the domain.}
\label{property_table}
\end{center}
\end{table}

The tests in \citet{guan2004isotropy} and \citet{maity2012test} require estimating the semivariogram and covariance function values, respectively, at the set of chosen lags. We denote the vector of point estimates of the semivariogram/covariance function at the chosen lags as $\widehat{\G}_n$, the true values as $\G$, and the asymptotic variance-covariance matrix of $\widehat{\G}_n$ as $\mathbf{\Sigma}$. Under increasing domain asymptotics, a central result for both methods is that
\begin{equation}\label{asymptoticMVN}
a_n(\widehat{\G}_n - \G) \xrightarrow[n \rightarrow \infty]{d} MVN(\mathbf{0}, \mathbf{\Sigma}),
\end{equation}
where $a_n$ is a normalizing constant. This result holds assuming stationarity of the RF and mild conditions on the RF's mixing and moments. See Section~\ref{discussion} for additional discussion regarding the assumptions for the tests. The test statistic is a quadratic form
\begin{equation}\label{ts}
TS = b_n^2(\A\widehat{\G}_n)^{\top}(\A\widehat{\mathbf{\Sigma}}\A^{\top})^{-1} (\A\widehat{\G}_n),
\end{equation}
where $\widehat{\mathbf{\Sigma}}$ is an estimate of the asymptotic variance-covariance matrix and $b_n$ is a normalizing constant. A $p$~value can be obtained from the asymptotic $\chi^2$ distribution with degrees of freedom given by the row rank of $\A$. Important differences between these works regard the distribution of sampling locations, shape of the sampling domain, and estimation of $\G$ and $\mathbf{\Sigma}$. Table~\ref{property_table} outlines the primary differences between these methods; we refer the interested reader to \citet{weller2015} for more details.



\section[Nonparametric tests implemented in spTest]{Nonparametric tests implemented in \pkg{spTest}}\label{sptest-functions}
The \R package \pkg{spTest} includes functions for implementing the tests developed in \citet{guan2004isotropy}, \citet{lu2005test}, and \citet{maity2012test}. The \pkg{spTest} functions for implementing these tests are listed in Table~\ref{property_table}. For example, the test from \citet{guan2004isotropy} for data observed at non-gridded, but uniformly distributed, sampling locations is implemented in the function \texttt{GuanTestUnif}, which takes the following arguments:

\begin{center}
\texttt{
GuanTestUnif(spdata, lagmat, A, df, h = 0.7, kernel = "norm",
  truncation = 1.5, xlims, ylims, grid.spacing = c(1, 1),
  window.dims = c(2, 2), subblock.h, sig.est.finite = T).
}
\end{center}

There are several necessary inputs. \texttt{spdata} must include the coordinates of sampling locations and the corresponding data values. This input can be a matrix or an object created by either the \pkg{sp} or \pkg{geoR} packages. The spatial lags used to estimate the semivariogram, denoted $\bL$, are specified in the matrix \texttt{lagmat}. The matrix $\A$ in \eqref{ts} is specified by \texttt{A} and provides the contrasts of the semivariogram estimates, and its row rank is indicated by the parameter \texttt{df} (the degrees of freedom for the asymptotic $\chi^2$ distribution). The values \texttt{h} and \texttt{kernel} provide the bandwidth (smoothing) parameter and form of the kernel smoother, respectively, used to smooth over spatial lags when estimating the semivariogram. If a normal smoothing kernel is used, then the \texttt{truncation} parameter indicates where to truncate the normal kernel (i.e., zero weight for spatial lags larger than this value). The parameters \texttt{xlims} and \texttt{ylims} give the horizontal and vertical limits of the sampling region (a rectangular sampling region is assumed). When performing a nonparametric tests of isotropy for non-gridded sampling locations, we must place a grid on the sampling domain and choosing a moving window or block size to estimate $\mathbf{\Sigma}$ in \eqref{asymptoticMVN} (see Section~\ref{nongrid-example}). A grid is placed over the sampling region according the width and height specified by \texttt{grid.spacing}. The dimensions of the moving window, given in the units of the underlying grid, are determined by the values in \texttt{window.dims}. The bandwidth of the smoothing kernel used to estimate the semivariogram on the subblocks of data created by the moving window is indicated by \texttt{subblock.h} and a finite sample adjustment to the estimate of the asymptotic variance-covariance matrix is made by setting \texttt{sig.est.finte = T}. For more information about the different arguments and guidelines on how to choose them, see \citet{weller2015}, the \pkg{spTest} manual, and the original works \citep{guan2004isotropy, guan2007asymptotic, maity2012test}.

\section[Applications: Using spTest to check for anisotropy]{Applications: Using \pkg{spTest} to check for anisotropy}\label{applications}
We demonstrate the functionality of the \pkg{spTest} on two data sets: the first containing data at gridded sampling locations; the second containing data collected via a non-gridded sampling design. For more details on the functions and examples using simulated data, see the \pkg{spTest} manual. The \pkg{spTest} package can be used independently of other packages built for analyzing spatial data, but it works nicely with two other packages loaded into \R: \pkg{fields} \citep{fields} and \pkg{geoR} \citep{geor}. We also load the \pkg{splines} \citep{rsoftware}, \pkg{MASS} \citep{mass}, and \pkg{rgdal} \citep{rgdal} packages, which we use to estimate mean functions, compute studentized residuals, and calculate map projection coordinates, respectively.


<<label = libs, echo = TRUE>>=
library("spTest")
library("geoR")
library("splines")
library("MASS")
library("rgdal")
library("mvtnorm")
@

For the two examples given below, we use graphical diagnostics and the hypothesis tests implemented in \pkg{spTest} to determine whether or not an assumption of isotropy is reasonable for spatially-referenced data. The general strategy will be to first do exploratory data analysis (EDA) of the original data and create a model for the mean of the spatial process using appropriate covariates. After estimating a model for the mean, we extract residuals and again use EDA to check for remaining spatial dependence and utilize graphical diagnostics and hypothesis tests to investigate potential anisotropy. For brevity, we have not included the full version of EDA code and plots; instead, we include only the most relevant to demonstrating the functionality of the \pkg{spTest} package. The complete version of the code is available on \texttt{github} \citep{wellerGIT}.

\subsection{Gridded sampling locations}\label{grid-section}
The gridded data used in this section come from the North American Regional Climate Change Assessment Program [NARCCAP] \citep{mearns2009regional}. The data set \texttt{WRFG} in \pkg{spTest} includes coordinates and a 24-year average of yearly average temperatures from runs of the Weather Research and Forecasting - Grell configuration (WRFG) regional climate model (RCM) using boundary conditions from the National Centers for Environmental Prediction (NCEP). The original data are available on the NARCCAP website and the \R code used to create the yearly averages is available on \texttt{github} \citep{wellerGIT2}. The data set contains both latitude and longitude and universal transverse mercator (UTM) coordinates. The UTM coordinates specify the regular grid for 14,606 grid boxes along with average temperature at surface at each grid box. Figure~\ref{narccap} displays a heat map of all of the data and was created using the \texttt{image.plot} function from the \pkg{fields} package . Due to computational considerations and because the methods in \pkg{spTest} assume stationarity, for our analysis we use a $20 \times 20$ subset of the grid boxes defined by the UTM coordinates over the central United States (see Figures~\ref{narccap} and ~\ref{hm1}).

<<label=narccap,fig=TRUE,include=FALSE>>=
data("WRFG")
coords <- expand.grid(WRFG$xc,WRFG$yc)
sub <- which(coords[,1] > 2.90e6 & coords[,1] < 3.95e6 
& coords[,2] > 1.2e6 & coords[,2] < 2.25e6)
coords.ll <- cbind((WRFG$lon-360)[sub], WRFG$lat[sub])
image.plot(WRFG$lon-360, WRFG$lat, WRFG$WRFG.NCEP.tas, 
col = two.colors(n = 256, start = "blue3", end = "red3", middle = "gray60"),
legend.lab = "Temp (K)",legend.cex = 0.8,legend.line = 2.2,
xlab = "Longitude", ylab = "Latitude", 
main = "Mean WRFG-NCEP Temperatures")
world(add = T)
left <- seq(1, 400, by = 20)
right <- seq(20, 400, by = 20)
for(i in 2:20){
	segments(coords.ll[i-1,1],coords.ll[i-1,2], coords.ll[i,1],coords.ll[i,2], lwd = 2)
	segments(coords.ll[left[i-1],1],coords.ll[left[i-1],2], coords.ll[left[i],1],coords.ll[left[i],2], lwd = 2)
	segments(coords.ll[right[i-1],1],coords.ll[right[i-1],2], coords.ll[right[i],1],coords.ll[right[i],2], lwd = 2)
	j <- i + 380
	segments(coords.ll[j-1,1],coords.ll[j-1,2], coords.ll[j,1],coords.ll[j,2], lwd = 2)
}
@

\begin{figure}
\begin{center}
\includegraphics[width = 3.75 in]{figures/fig-narccap.eps}
\caption{Heat map of temperatures in latitude and longitude coordinates. The temperature values are the 24 year average of WRFG-NCEP yearly average temperature. The black box indicates the relevant subset of data used for the example in Section~\ref{grid-section}, and this subset is displayed in Figure~\ref{hm1}. }
\label{narccap}
\end{center}
\end{figure}

To investigate potential anisotropy in the relevant subset of these data, we can examine two graphical diagnostics: a heat map and directional sample semivariograms. We use the function \texttt{variog4} from the \pkg{geoR} package to estimate directional semivariograms to visually assess isotropy properties.

<<label = heatmap1, results = hide, fig = T, include = F, width = 6>>=
tas <- c(WRFG$WRFG.NCEP.tas)[sub]
x.coord <- unique(coords[sub,1])
y.coord <- unique(coords[sub,2])
nx <- length(x.coord)
ny <- length(y.coord)
tas.mat <- matrix(tas, nrow = nx, ncol = ny, byrow = F)
image.plot(x.coord, y.coord, tas.mat, 
col = two.colors(n = 256, start = "blue3", end = "red3", middle = "gray60"),
legend.lab = "Temp (K)",legend.cex = 0.8,legend.line = 2.2,
ylab = "Northing", xlab = "Easting",
main = "Subset of Temperatures")
@
<<label = dirsemi1, results = hide, fig = T, include = F>>=
tas.geodat <- as.geodata(cbind(coords[sub,1],coords[sub,2], tas))
plot( variog4(tas.geodat), xlab = "distance (meters)", ylab = "estimated semivariogram")
title("Directional Sample Semivariograms")
@
\begin{figure}
	\begin{subfigure}[t]{0.5\linewidth}
	\centering
	\includegraphics[width = 2.5 in, height = 2.5 in]{figures/fig-heatmap1.eps}
	\caption{}
	\label{hm1}
	\end{subfigure}
	\begin{subfigure}[t]{0.5\linewidth}
	%\centering
	\includegraphics[width = 2.5 in]{figures/fig-dirsemi1.eps}
	\caption{}
	\label{dv1}
	\end{subfigure}
\caption{Graphical assessments of isotropy in the $20 \times 20$ subset of WRFG temperature data. Because northing coordinates have not been accounted for, the heat map (Figure~\ref{hm1}) indicates that the dependence between observations is stronger in the east-west direction than the north-south direction. The directional dependence is also evidenced by the differences between the directional sample semivariograms (Figure~\ref{dv1}).}
\label{fig2}
\end{figure}

The heat map in Figure~\ref{hm1} indicates that the spatial process is anisotropic, having a stronger spatial dependence in the horizontal direction than the vertical direction. Intuitively, northing coordinates are an important factor in determining average temperature, and we need to include its effect in a model for these data. We also notice non-linear trends in temperature as a function of easting coordinates in Figure~\ref{hm1}. Thus, the anisotropy can be attributed, at least in part, to the fact that we have not modeled important covariates related to the process. The directional sample semivariograms in Figure~\ref{dv1} reaffirm the notion that the data exhibit anisotropy as the $90^{\circ}$ sample semivariogram appears much different than the other three. Before modeling the effects of northing and easting coordinates, we use the \texttt{GuanTestGrid} function from \pkg{spTest} to affirm our understanding that these data exhibit anisotropy.

Necessary conditions for the asymptotic properties of the nonparametric tests to hold are typically met when the data are Gaussian (see Section~\ref{discussion}). A quantile-quantile (QQ) plot (not shown) of the relevant subset of WRFG temperatures indicates that a Gaussian assumption is reasonable. To implement the nonparametric test in \citet{guan2004isotropy} via the function \texttt{GuanTestGrid}, we need to specify the spatial lags that will be used to test for differences in the semivariogram. For this test we choose the lag set \eqref{lagset} and use the matrix $\A$ in \eqref{amat} to contrast the semivariogram estimates. With the first row of $\A$ and the first two entries of $\bL$, we are contrasting the estimated dependence structure in the $0^{\circ} (\h_1)$ and $90^{\circ} (\h_2)$ directions for data separated by one horizontal or vertical sampling location. The second row of $\A$ and second two entries of $\bL$ contrast the estimated dependence structure in the $45^{\circ} (\h_3)$ and $135^{\circ} (\h_4)$ directions for data separated by one diagonal sampling location. Because the grid spacing between sampling locations is 50,000 meters, we set the the scaling parameter \texttt{delta = 50,000}. To create subblocks of data used to estimate $\mathbf{\Sigma}$ in \eqref{ts}, we choose a moving window with a size of $4 \times 4$ grid cells. The moving window dimensions should be chosen so that the window has the same shape (i.e., square or rectangle) and orientation as the sampling domain. To maximize the amount of data used to estimate $\mathbf{\Sigma}$, the dimensions of the window should evenly divide the number of columns and rows, respectively, of the entire region. The window dimensions should also be compatible with the spatial lags in $\bL$. For example, if sampling locations are on the integer grid $\mathbb{Z}^2$, a window with dimensions of $2 \times 2$ grid cells cannot be used to estimate the variability of the semivariogram at a lag with Euclidean distance longer than $\sqrt{2}$, the maximum distance between locations in the moving window. For this example there are 20 rows and columns, and we are using lags with spacings of one or two sampling locations; hence, window dimensions of $2 \times 2$ or $4 \times 4$ grid cells are reasonable choices. We run the test using window dimensions of $4 \times 4$ grid cells via the following code, suppressing some of the output for brevity.

<<label = test1 >>=
my.delta <- 50000
mylags <-  rbind(c(1,0), c(0, 1), c(1, 1), c(-1,1))
myA <-  rbind(c(1, -1, 0 , 0), c(0, 0, 1, -1))
tr <- GuanTestGrid(spdata = tas.geodat, delta = my.delta, lagmat = mylags, A = myA,
 df = 2, window.dims = c(4,4), pt.est.edge = TRUE, sig.est.edge = TRUE, 
 sig.est.finite = TRUE)
tr$alternative <- NULL
tr$sigma.hat <- NULL
print(tr)
@

As we suspected, the results of the hypothesis test ($p$~value $< 0.05$) indicate that the data exhibit anisotropy. We note that for gridded data, \citet{guan2004isotropy} recommend using the $p$~value computed via a finite sample correction. The function \texttt{GuanTestGrid}, and other functions in \pkg{spTest}, return a $p$~value(s) for the test and information used in computing the $p$~value, such as the point estimates ($\widehat{\G}_n$), estimates of the asymptotic variance-covariance matrix ($\widehat{\mathbf{\Sigma}}$), the number of subblocks used to estimate $\mathbf{\Sigma}$, and other information about the estimation process. We note that the point estimates for the directional semivariograms are slightly different between the functions from the \pkg{spTest} and \pkg{geoR} packages due to different kernel methods used in estimation.

As previously mentioned, we need to model the effects of northing and easting UTM coordinates on average temperature. We fit temperature as a nonparametric additive function of both the northing and easting coordinates via least-squares using cubic splines. The cubic splines can be specified using the function \texttt{ns} from the \pkg{splines} package and the least squares fit is computed via the \texttt{lm} function.

<<label = model1, results = hide>>=
m1 <- lm(tas ~ ns(coords[sub,1], df = 3) + ns(coords[sub,2], df = 3))
summary(m1)
@

After removing the mean effects of the coordinates, we can check for any remaining (unaccounted for) spatial dependence and evidence of anisotropy in the residuals using graphical diagnostics and a hypothesis test. A QQ plot of the studentized residuals (not shown) indicates that a Gaussian assumption is reasonable.

<<label = iso2, results = hide, fig = T, include = F, width = 6>>=
resid.mat <- matrix(studres(m1), nrow = nx, ncol = ny, byrow = F)
image.plot(x.coord, y.coord, resid.mat, 
col = two.colors(n = 256, start = "blue3", end = "red3", middle = "gray60"),
xlab = "Easting", ylab = "Northing")
title("Heat Map of Studentized Residuals")
@
<<label = dsvar2, results = hide, fig = T, include = F>>=
resid.geo <- as.geodata(as.matrix(cbind(coords[sub,1:2], studres(m1))))
plot(variog4(resid.geo), xlab = "distance (meters)", ylab = "estimated semivariogram")
title("Directional Sample Semivariograms")
@
\begin{figure}
	\begin{subfigure}[t]{0.5\linewidth}
	\centering
	\includegraphics[width = 2.5 in, height = 2.5 in]{figures/fig-iso2.eps}
	\caption{}
	\label{hm2}
	\end{subfigure}
	\begin{subfigure}[t]{0.5\linewidth}
	%\centering
	\includegraphics[width = 2.5 in]{figures/fig-dsvar2.eps}
	\caption{}
	\label{dv2}
	\end{subfigure}
\caption{Graphical assessments of isotropy in the studentized residuals from the WRFG temperature data. The appearance of elongated areas of similar residual values in the heat map (Figure~\ref{hm2}) indicates that the process may be anisotropic. The directional semivariograms (Figure~\ref{dv2}) do not appear to exhibit differences, indicating that the process is isotropic. A nonparametric test of isotropy can assist in determining whether or not an assumption of isotropy is reasonable.}
\label{fig3}
\end{figure}

The clusters of similar values in the heat map of Figure~\ref{hm2}, and the increase, followed by a leveling off, of the semivariogram values as distance increases in the directional sample semivariograms in Figure~\ref{dv2} indicate that the residuals are still spatially dependent. However, the plots in Figure~\ref{fig3} do not clearly illustrate whether or not the residuals exhibit anisotropy. There appears to be directional dependence along the NW to SE direction in the northern parts of the heatmap (Figure~\ref{hm2}). The directional sample semivariograms do not appear to be different until the distance is greater than 200,000 meters. Semivariogram estimates at large distances can be unreliable and a general recommendation is to compute the semivariogram for lags with distance that is less than half of the maximum observed distance separating sampling locations \citep[pg.~155]{schab}. Additionally, directional semivariograms are less reliable than a uni-directional semivariogram because fewer pairs of sampling locations are used at each distance for directional estimation. The unreliability of the sample semivariograms at the larger distances, coupled with the lack of a measure of uncertainty, make it difficult to determine whether or not an assumption of isotropy is reasonable using a plot of the directional sample semivariogams. To gain more insight into the isotropy properties, we perform a nonparametric hypothesis test of isotropy using the residuals with the same choices for $\bL$, $\A$, and the window dimensions.

<<label = test2>>=
tr <- GuanTestGrid(spdata = resid.geo, delta = my.delta, lagmat = mylags, A = myA, df = 2, window.dims = c(4,4))
tr$p.value.finite
@

Here the residuals do not provide evidence for anisotropy ($p$~value $> 0.05$). These results suggest that it may be appropriate to choose an isotropic covariance function to model the residuals. However, it is important to note that we have not included the effect of other potentially influential covariates such as elevation or water cover in the model for temperature. Additionally, although we examined a $20 \times 20$ subset of the data, the grid boxes still cover a large geographic region of the U.S., and thus an assumption of stationarity, which is needed for the asymptotic properties of the hypothesis test to hold, may not be reasonable.

\subsection{Non-gridded sampling locations}\label{nongrid-example}
The non-gridded data set used in this section describes monthly surface meterology in a region of the state of Colorado and comes from the National Center for Atmospheric Research (NCAR). The data are available in the \R package \pkg{fields}. For this example, our variable of interest is the log (mm) of the 30-year average of average yearly precipitation at 344 station locations during the time period 1968-1997. 

Like the temperature data, our goal will be to model the mean effect of covariates and check for spatial dependence and potential anisotropy in the residuals. The first two steps of the analysis are to compute the yearly precipitation averages and convert the latitude/longitude coordinates to UTM coordinates. We divide the UTM coordinates by 100,000 so that distances are measured in hundreds of kilometers. Scaling the coordinates eases the choice of tuning parameters for the test in \citet{guan2004isotropy}. To meet the Gaussian assumption, we take the log transform of the average precipitation. Figure~\ref{qplots} displays quilt plots of the log precipitation values and the elevation of the stations.

<<label = rain, results = hide, fig = T, include = F>>==
data("COmonthlyMet")
sub30 <- CO.ppt[74:103,,]
nstations <- 376
years <- 1968:1997
nyears <- length(years)
yr.avg <- matrix(data = NA, nrow = nstations, ncol = nyears)
for(i in 1:nyears){
	yr.dat <- sub30[i,,]
	yr.avg[,i] <- apply(yr.dat, 2 , mean, na.rm = T)
}
avg30 <- apply(yr.avg, 1, mean, na.rm = T)
CO.loc <- as.matrix(CO.loc)
CO.loc.utm <- project(CO.loc, 
"+proj=utm +zone=13 ellps=WGS84")/1e+5
quilt.plot(CO.loc.utm, log(avg30), 
col = two.colors(n = 256, start = "blue3", end = "red3", middle = "gray60"),
legend.lab = "Precip (log mm)",legend.cex = 0.8,legend.line = 2.2,
xlab = "Easting", ylab = "Northing",
main = "Quilt Plot of log(precip)")
mp <- map("state", region = c("colorado", "wyoming","nebraska", "utah", "new mexico", "oklahoma"), plot = F)
st.loc <- cbind(mp$x, mp$y)
states <- project(st.loc, 
"+proj=utm +zone=13 ellps=WGS84")/1e+5
points(states[,1], states[,2], type = "l", lwd = 1.5)
@
<<label = elev, results = hide, fig = T, include = F>>==
quilt.plot(CO.loc.utm, CO.elev, 
col = two.colors(n = 256, start = "blue3", end = "red3", middle = "gray60"),
legend.lab = "Elevation (meters)",legend.cex = 0.8,legend.line = 2.7,
xlab = "Easting", ylab = "Northing",
main = "Quilt Plot of Elevation")
points(states[,1], states[,2], type = "l", lwd = 1.5)
@

\begin{figure}
	\begin{subfigure}[t]{0.5\linewidth}
	\centering
	\includegraphics[width = 2.75 in]{figures/fig-rain.eps}
	\caption{}
	\label{qp-precip}
	\end{subfigure}
	\begin{subfigure}[t]{0.5\linewidth}
	%\centering
	\includegraphics[width = 2.75 in]{figures/fig-elev.eps}
	\caption{}
	\label{qp-elev}
	\end{subfigure}
\caption{Quilt plots showing the locations of the weather stations in Colorado and surrounding region along with the log of average yearly precipitation (\ref{qp-precip}) and elevation (\ref{qp-elev}) at each station.}
\label{qplots}
\end{figure}

Colorado has two distinct geographic regions: the mountainous region in the west and the plains region in the east. Figure~\ref{qp-elev} illustrates these two regions, and we can begin to notice a possible relationship between elevation and average precipitation. We explore the potential relationship between log precipitation and elevation using scatter plots (see Figure~\ref{elev_scatter}).

<<label = rain1, results = hide, fig = T, include = F>>==
plot(CO.elev, log(avg30) , xlab = "Elevation (meters)",ylab = "Precip (log mm)",main = "log(Precip) vs. Elevation")
m1 <- lm( log(avg30) ~ ns(CO.elev, df = 3))
summary(m1)
fits <- m1$fitted.values
bad <- is.na(avg30)
x <- CO.elev[which(!bad)]
lines(sort(x), fits[order(x)], lwd = 3, col = "red")
@
<<label = rain2, results = hide, fig = T, include = F>>==
qqnorm(studres(m1))
abline(0,1)
@

\begin{figure}
	\begin{subfigure}[t]{0.5\linewidth}
	\centering
	\includegraphics[width = 2.5 in]{figures/fig-rain1.eps}
	\caption{}
	\label{elev_scatter}
	\end{subfigure}
	\begin{subfigure}[t]{0.5\linewidth}
	%\centering
	\includegraphics[width = 2.5 in]{figures/fig-rain2.eps}
	\caption{}
	\label{rqq}
	\end{subfigure}
\caption{Results from the model relating log(precipitation) and elevation. Figure~\ref{elev_scatter} displays the nonparametric fit relating elevation to log(precipitation). Figure~\ref{rqq} shows the QQ plot of studentized residuals from the nonparametric fit.}
\label{rain_model}
\end{figure}

We fit a cubic smoothing spline via least squares to model the relationship between log(precipitation) and elevation. The estimate is shown in Figure~\ref{elev_scatter}, and a QQ plot of residuals in Figure~\ref{rqq} indicates that a Gaussian assumption is reasonable. We will use the residuals from this model to check for remaining spatial dependence and potential anisotropy. We use \texttt{variog4} to estimate directional sample semivariograms. 

<< label = precipdv, results = hide, fig = T, include = F>>=
precip.resid <- cbind(CO.loc.utm[which(!bad),][,1],CO.loc.utm[which(!bad),][,2] , studres(m1))
precip.geo <- as.geodata(precip.resid) 
plot(variog4(precip.geo), xlab = "distance (100s of km)", ylab = "estimated semivariogram", legend = F)
legend("bottomright", legend = c(expression(0*degree), expression(45*degree), expression(90*degree), expression(135*degree)), col = 1:4, lty = 1:4)
title("Directional Sample Semivariograms")
@

\begin{figure}[h]
\begin{center}
\includegraphics[width = 3.5 in]{figures/fig-precipdv.eps}
\end{center}
\caption{A graphical assessment of isotropy in the studentized residuals from the model relating log(precipitation) to elevation generated by using the \pkg{geoR} function \texttt{variog4}. The directional sample semivariogram in the $0^{\circ}$ direction appears to be different from the other three at distances less than 2. Because there is no measure of uncertainty, it can be difficult to determine whether or not an assumption of isotropy is reasonable.}
\label{raindv}
\end{figure}

The increase, followed by a leveling off, of the semivariogram values as distance increases in Figure~\ref{raindv} indicates that there is spatial dependence remaining in the data. We notice that the $0^{\circ}$ semivariogram appears to be slightly different than the other three, but there is no measure of uncertainty, so we cannot determine if the differences are statistically significant. The sample semivariograms also suggest the presence of a large amount of small scale variation, often called a nugget effect, in the data. This can be seen by noting the ratio of semivariogram values at the shortest observed distance to the semivariogram values at the longest reliable distance. The large jumps and decrease in the estimated semivariogram values in Figure~\ref{raindv} indicate that semivariogram estimates become unreliable beyond a distance of two. The sample semivariogram values in Figure~\ref{raindv} are approximately 0.7 at the shortest distance and approximately 1.1 at a distance of two. This suggests that approximately $63\% \approx (0.7/1.1)100\%$ of the variability in the data is due to small scale variation. We note the apparently large nugget effect because this small scale variation is detrimental to the size and power of nonparametric tests of isotropy \citep{weller2015}. Despite the small scale variation, we will proceed with nonparametric hypothesis tests to assist in determining if an assumption of isotropy is reasonable. 

There are two procedures for testing isotropy in non-gridded data available in \pkg{spTest}: \citet{guan2004isotropy} and \citet{maity2012test}. To choose between these two, we need to decide whether or not it is reasonable to assume that sampling locations are uniformly distributed on the sampling domain. The methods for non-gridded data from \citet{guan2004isotropy} rely on the assumption that sampling locations are uniformly distributed while \citet{maity2012test} can be used on any general sampling design. To check this assumption, we can turn to methods from the spatial point process literature to perform a test of complete spatial randomness (CSR) (i.e., a uniform spatial distribution) for the sampling locations. Methods for testing CSR are available in the \R package \pkg{spatstat} \citep{baddeley2004spatstat}. For brevity, we do not display the results of the CSR test here, but note that they do not provide evidence against the assumption of CSR for these sampling locations so either test of isotropy can be used.

For both \citet{guan2004isotropy} and \citet{maity2012test}, we need to choose the lag set, $\bL$, and the contrast matrix, $\A$. Because semivariogram estimates appear to be unreliable at distances greater than two, we should choose lags having Euclidean distance less than this distance. We choose the lag set 
\[
 \bL = \{ \h_1 = (0.60,0), \h_2 = (0,0.60), \h_3 = (0.45,0.45),  \h_4 = (-0.45, 0.45)\},
\]
and again we use the matrix $\A$ in \eqref{amat}. The set $\bL$ corresponds to lags in the $0^{\circ}, 90^{\circ}, 45^{\circ}, 135^{\circ}$ directions, respectively, having Euclidean distances of $||\h_1|| = ||\h_2|| = 0.60$ (60 km) and $||\h_3|| = ||\h_4|| \approx 0.64$ (64 km).
<<label = lags, keep.source = TRUE>>=
mylags <- rbind(c(0.60, 0), c(0, 0.60), c(0.45, 0.45), c(-0.45, 0.45))
myA <- rbind(c(1, -1, 0 , 0), c(0, 0, 1, -1))
@

The next step in implementing the methods from \citet{guan2004isotropy} and \citet{maity2012test} is to determine the size of the moving windows and the block size, respectively, used to estimate the asymptotic variance-covariance matrix, $\mathbf{\Sigma}$. The moving window is shifted over the sampling region, creating subblocks of data used to estimate $\mathbf{\Sigma}$. Likewise, for the test in \citet{maity2012test}, the block size is used to implement the grid-based block bootstrap [GBBB] \citep{lahiri2006resampling}.

There are two steps in determining the appropriate window/block size for non-gridded sampling locations. First, we place a grid over the sampling domain; second, we specify scaling parameters that will define the window/block size in terms of that grid. We should complete this two step process while keeping three goals in mind: (1) the number of sampling locations per window/block, denoted $n_b$, should be approximately $\sqrt{n}$ \citep{weller2015}; (2) the windows/blocks should have have the same orientation (i.e., square or rectangle) as the entire sampling domain; and (3) the scaling parameters should be compatible with the dimensions of the underlying grid. 

For the Colorado precipitation data, recall that one unit of distance equals 100 km. The dimensions of the sampling region are approximately $7.3 \times 5.5$ (width $\times$ height), providing a total area of 40.15. For $n = 344$ uniformly distributed sampling locations, we expect approximately $344/40.15 = 8.6$ sampling locations per unit area. Recalling goal (1), we seek to construct windows/blocks with $n_b \approx \sqrt{344} = 18.5$ sampling locations, or equivalently, windows/blocks with an area of approximately $18.5/8.6 = 2.15$. Goal (2) indicates we want to create rectangular windows/blocks with slightly larger width than height, and (3) says that if our grid divides the $x$-axis into 12 grid cells, then the scaling parameter defining the width of the window/block should be 3 or 4 because those numbers evenly divide 12. For the CO precipitation data, if we choose our grid to divide the $x$-axis into 16 cells and the $y$-axis into 12 cells, we have a grid with $(x,y)$ spacing of roughly $(7.3/16, 5.5/12) \approx (0.46, 0.46)$. The resulting grid is plotted in Figure~\ref{gridplot}. Then, choosing our scaling parameters to be $4 \times 3$, we have windows/blocks with dimensions of approximately $(4)(0.46) \times (3)(0.46) = 1.84 \times 1.38$ and area of $(1.84)(1.38) \approx 2.54$, or equivalently with an expected number of points per block of $n_b = (2.54)(8.6) \approx 21.8$.

<<label = grid, fig = T, include = F>>==
my.color <- two.colors(n = 256, start = "blue3", end = "red3", middle = "gray60")
quilt.plot(precip.resid[,1:2], precip.resid[,3], 
col = my.color,
xlab = "Easting", ylab = "Northing", xlim = c(0.75,8.65), ylim = c(40.1, 46.2) )
title("Quilt Plot of Residuals and Grid Used for Subsampling")
tol <- 0.02
my.xlims <-  c(min(precip.resid[,1]) - tol, max(precip.resid[,1]) + tol )
my.ylims <-  c( min(precip.resid[,2]) - tol, max(precip.resid[,2]) + tol )
xlen <-  my.xlims[2]-my.xlims[1]
ylen <-  my.ylims[2]-my.ylims[1]
my.grid.spacing <-  c(xlen/16, ylen/12)
xgrid <- seq(my.xlims[1], my.xlims[2], by = my.grid.spacing[1])
ygrid <- seq(my.ylims[1], my.ylims[2], by = my.grid.spacing[2])
segments(x0 = xgrid, y0 = min(my.ylims), y1 = max(my.ylims), lty = 2)
segments(x0 = min(my.xlims), y0 = ygrid, x1 = max(my.xlims), lty = 2)
@

\begin{figure}[h]
\begin{center}
\includegraphics[width = 3.55 in]{figures/fig-grid.eps}
\end{center}
\caption{Quilt plot of the studentized residuals from the model relating elevation to log(precipitation) for the weather station locations. The grid placed on the region that is used to define the moving windows in \citet{guan2004isotropy} and block size in \citet{maity2012test} is also shown. Because the sampling locations are not gridded, it can be difficult to assess isotropy properties via a quilt plot.}
\label{gridplot}
\end{figure}

For the functions \texttt{GuanTestUnif} and \texttt{MaityTest}, the upper and lower limits of the sampling region in the $x$ and $y$ directions are given by the \texttt{xlims} and \texttt{ylims} arguments. Note that the values defining the upper and lower limits should be slightly larger than the minimum and maximum observed $x$ and $y$ coordinates. The horizontal and vertical spacing, respectively, of the grid laid on the sampling region is defined by the two values in \texttt{grid.spacing}. The horizontal and vertical scaling parameters that define the size of the moving windows in \texttt{GuanTestUnif} and blocks in \texttt{MaityTest} in terms of the underlying grid are given by the \texttt{window.dims} and \texttt{block.dims} arguments, respectively. We recommend using visualizations of different grid choices and algebraic calculations, as done above, to assist in choosing a grid and the window/block dimensions. When the scaling parameters defining the moving window or block dimensions are not compatible with the number of rows or columns of gridded sampling locations or the dimensions of the grid laid on the sampling region for non-gridded locations, the functions in \pkg{spTest} will print a warning message because they do not currently handle partial (incomplete) blocks. Likewise, if the chosen window or block dimensions for non-gridded sampling locations creates (sub)blocks of data with few or no sampling locations, the functions \texttt{GuanTestUnif} and \texttt{MaityTest} will discard (sub)blocks that do not have enough sampling locations and print a warning message. The $p$~value of the hypothesis test will be sensitive to the choice of moving window and block dimensions. See \citet{weller2015} and the original works \citep{guan2004isotropy, maity2012test} for more recommendations on choosing these values.

The next step for implementing the test in \citet{guan2004isotropy} is choosing the smoothing (bandwidth) parameters for smoothing over lags on the entire domain and within each subblock created by the moving windows. The smoothing parameters should be chosen based on the number and density of the sampling locations with larger values of the smoothing parameter inducing higher levels of smoothing, i.e., allowing a greater distance and direction tolerance. In our experience, smoothing parameter values between 0.6 and 0.9 tend to produce reasonable results when using a standard normal Gaussian smoothing kernel truncated at 1.5. However, the $p$~value of the hypothesis test will change with the bandwidth. For this example, we choose a bandwidth of \texttt{h = 0.70} for smoothing over lags on the entire domain, and a bandwidth of \texttt{subblock.h = 0.85} for smoothing over lags on on the subblocks of data created by the moving window. Choosing a larger bandwidth for the subblocks equates to allowing for a larger lag distance and direction tolerance, which is needed for subblocks that have few sampling locations. We also use the default Gaussian smoothing kernel (\texttt{kernel = "norm"}) truncated at 1.5 (\texttt{truncation = 1.5}). Because the sample size is less than 500, we use a finite sample adjustment to approximate the $p$ value \citep{guan2004isotropy, weller2015}.

Finally, for the test in \citet{maity2012test} we need to choose the number of bootstrap resamples that will be used in the GBBB procedure to estimate $\mathbf{\Sigma}$. We recommend using at least 50 bootstrap samples; however, the bootstrapping procedure is computationally intensive. We choose \texttt{nBoot = 100} bootstrap samples for our example, and we note that the number of bootstraps does not affect the precision of the $p$~value, which is computed via the asymptotic $\chi^2$ distribution. Having determined values for the different options, we can now perform the hypothesis tests. For reproducibility of the bootstrap in the \texttt{MaityTest} function, we set the random seed.

<<label = nongrid>>=
myh <-  0.7
myh.sb <-  0.85
tr.guan <- GuanTestUnif(spdata = precip.resid, lagmat = mylags, A = myA, 
df = 2, h = myh, kernel = "norm", truncation = 1.5, xlims = my.xlims, 
ylims = my.ylims, grid.spacing = my.grid.spacing, window.dims = c(4,3), 
subblock.h = myh.sb)
tr.guan$p.value.finite

set.seed(2016)
tr.maity <- MaityTest(spdata = precip.resid, lagmat = mylags, 
A = myA, df = 2, xlims = my.xlims, ylims = my.ylims,
grid.spacing = my.grid.spacing, block.dims = c(4,3), nBoot = 100)
tr.maity$p.value
@

For both of the tests, the data provide evidence in favor of anisotropy ($p$~value $<$ 0.05). Thus, an isotropic model may be appropriate for modeling the residuals. Additionally, the apparent anisotropy may also be present due to unaccounted for covariates (e.g., northing/easting coordinates).

\section{Discussion}\label{discussion}
Choosing a covariance function is an important step in modeling spatially-referenced data and a variety of choices for the covariance function are available (e.g., anisotropy, nonstationarity, parametric forms). The \R package \pkg{spTest} implements several nonparametric tests for checking isotropy properties which avoid specifying a parametric form for the covariance function. \citet{weller2015} perform a simulation study comparing the empirical size and power of the methods for different degrees of anisotropy. They find that methods from \citet{guan2004isotropy} tend to outperform the competitor for gridded and non-gridded data.

One concern regarding the methods in \pkg{spTest} is that they tend to have low power when the anisotropy is weak and the data are not gridded \citep{weller2015, guan2004isotropy, maity2012test}. A second concern is that the results of the tests are potentially sensitive to user choices, for example, the moving window size and bandwidth in the method from \citet{guan2004isotropy}. The optimal choices for these values is still an open question. \citet{weller2015} offer further recommendations for how to choose the user defined values, such as the window size and bandwidth, based on simulated data. Finally, as noted earlier, the size and power of the methods are adversely affected by the presence of small scale variation (nugget effect). Because of these concerns, we recommend using the nonparametric methods in conjunction with graphical techniques.

An implicit assumption of the methods discussed in this paper is ergodicity of the spatial process, an assumption that is difficult to verify \citep[pg.~57-58]{cressie1993}. However, there are two important assumptions which practitioners should consider. The first is an assumption of strict stationarity, \eqref{strict}. While this assumption is difficult to check, it follows from assuming the RF is weakly stationary and Gaussian \eqref{stationary_rel}. The assumption of Gaussian data lies at the heart of many spatial analyses \citep{gelfand2016spatial} and is easily checked with a QQ plot. The assumption of weak stationarity may be questionable for spatial data over large geographic regions and methods have been developed for testing this assumption \citep[see, e.g., ][]{corstanjeNS, fuentesNStest, jun2012test, bandyopadhyay2015test}. The second assumption required is a mixing condition that states the dependence between observations goes to 0 at large distances \citep{hall1994properties, sherman1994nonparametric}. In the case of a stationary GRF, this condition is met when the covariance function, $C(\h)$, is 0 for sufficiently large $||\h||$, which also implies ergodicity \citep[pg.~58]{cressie1993}. One way to check this assumption in practice is to look for a leveling off of the sample covariogram or semivariogram values as distance increases, indicating that the are data nearly independent at large distances.

After determining whether or not an assumption of isotropy is reasonable, we can choose a parametric or nonparametric model for the covariance function. \citet{weller2015} further illustrate the role of nonparametric tests of isotropy in the process of modeling spatially-referenced data. We have demonstrated how graphical techniques and the functions available in the \R package \pkg{spTest} can be used in a complementary role to check for anisotropy. Future work includes extending the functionality of \pkg{spTest} to handle non-rectangular sampling domains and improving computational efficiency by programming functions in \proglang{C++}.

\section*{Acknowledgements}
Weller's work was partially supported by the National Science Foundation Research Network on Statistics in the Atmospheric and Ocean Sciences (STATMOS) through grant DMS-1106862 and AGS-1419558. We wish to thank the North American Regional Climate Change Assessment Program (NARCCAP) for providing the data used in this paper. NARCCAP is funded by the National Science Foundation (NSF), the U.S. Department of Energy (DoE), the National Oceanic and Atmospheric Administration (NOAA), and the U.S. Environmental Protection Agency Office of Research and Development (EPA). The author would like to thank members of the Hooten-Hoeting lab and three anonymous reviewers for helpful comments that improved the manuscript.

\newpage

\bibliographystyle{apa}
\bibliography{spTestReferences}

\end{document}
