# Robust Maximum L$q$-Likelihood Covariance Estimation
for Replicated Spatial Data

Sihan Chen  Joydeep Chowdhury  and Marc G. Genton
Statistics Program, King Abdullah University of Science and Technology
{sihan.chen, joydeep.chowdhury<EMAIL_ADDRESS>
###### Abstract

The maximum Lq-likelihood estimator (MLqE) is a parameter estimator
alternative to the maximum likelihood estimator (MLE) that considers the
$q$-th power of the likelihood values for some $q<1$. In this method, extreme
values are down-weighted because of their lower likelihood values, which
yields robust estimates. In this work, we study the properties of the MLqE for
spatial data with replicates. We investigate the asymptotic properties of the
MLqE for Gaussian random fields with a Matérn covariance function, and carry
out simulation studies to investigate the numerical performance of the MLqE.
We show that it can provide more robust and stable estimation results when
some of the replicates in the spatial data contain outliers. In addition, we
develop a mechanism to find the optimal choice of the hyper-parameter $q$ for
the MLqE. The robustness of our approach is further verified on a United
States precipitation data set. Compared with other robust methods for spatial
data, our proposal is more intuitive and easier to understand, yet it performs
well when dealing with data sets containing outliers.

##  1 Introduction

The Maximum Likelihood Estimator (MLE, aldrich1997ra) has been one of the most
powerful statistical method widely used in various domains. Its applications
expanded across various fields, including econometrics (greene1980maximum;
cramer1989econometric), genetics (shaw1987maximum; beerli2006comparison) and
geography (elhorst2005unconditional; ree2008maximum), and gained popularity
for its asymptotic properties, such as consistency and efficiency (see, e.g.,
casella2021statistical). MLE also plays an important role while analyzing
spatial and spatio-temporal data (mardia1984maximum; moller2008handbook), and
in recent years, there are some tools developed for efficiently computing the
MLE for such data, for example ExaGeoStat (exageostat) and GpGp
(katzfuss2021general).

However, while using the MLE with datasets containing outliers, the likelihood
value might be largely distorted by these extreme values, and this may lead to
inaccurate estimation results (chen2014comparison). The maximum
L$q$-Likelihood estimator (ML$q$E) proposed by lq provides an effective way to
tackle this problem. Unlike the MLE which calculates the direct summation of
the likelihood values at each data point, the ML$q$E sums the $q$-th power of
the likelihood values with some $q<1$, which down-weights the extreme values
and leads to robust estimation. lq and ferrari established the asymptotic
normality of the ML$q$E under some regularity conditions and provided the
expression for its asymptotic variance.

In this work, we aim to employ the ML$q$E method for spatio-temporal data.
Here, we consider the spatial data at different time points as realizations of
a random field across the spatial domain. We focus on Gaussian random fields
with Matérn covariance function, which is one of the most influential model
for spatial statistics. The Matérn covariance function was proposed by
matern1960spatial, and later became popular with handcock1993bayesian. It is
applied in many of the most well known R packages, such as RandomFields
(schlather2013randomfields) and fields (nychka2021fields). A comprehensive
review of the application of the Matérn covariance function in spatial
statistics and many other related research fields can be found in
porcu2023mat. There are several different widely accepted ways to parameterise
the Matérn covariance function; see wang2023parameterization for a
comprehensive comparison of the three most popular parameterisations. In this
article, we adopt the following formulation of the Matérn covariance function,
which is the first type of parametrisation in wang2023parameterization:

$\displaystyle\mathcal{M}(h;\bm{\theta})=\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left(\frac{h}{\beta}\right)^{\nu}K_{\nu}\left(\frac{h}{\beta}\right).$
(1.1)

Here, $h$ is the distance between a pair of spatial locations, and
$\bm{\theta}=(\sigma^{2},\beta,\nu)^{\top}$ is the parameter vector that we
aim to estimate in a robust way with our proposed method. Its components
$\sigma^{2},\beta,\nu$ are the variance, range, and smoothness parameters
respectively, and $\mathcal{K}_{\nu}(\cdot)$ is the modified Bessel function
of the second kind of order $\nu$.

Figure 1: In the first sub-figure we present the functional boxplot of the
variograms of the January precipitation data, which we use to identify the
outliers in the dataset. The remaining three sub-figures present MLE and
ML$q$E estimation results of the three parameters, variance
($\hat{\sigma}^{2}$), range ($\hat{\beta}$) and smoothness ($\hat{\nu}$), with
the January data from 1928 to 1997. The black curves represent the ML$q$E
results with different values of $q$, and when $q=1$ (indicated by the blue
dot on the curve), the ML$q$E is the same as the MLE. The blue horizontal
lines and the red dashed lines denote the MLE results with and without the
outliers in the dataset.

We first demonstrate how the ML$q$E diminishes the influence of the outliers
in the dataset on the parameter estimation results using a motivating example
with the US precipitation data retrieved from
https://www.image.ucar.edu/Data/US.monthly.met/USmonthlyMet.shtml. The
detailed definition of the ML$q$E can be found in Section 2. The dataset that
we consider here contains hourly precipitation data from 1895 to 1997 recorded
at thousands of observation stations throughout the country. In this example,
we consider the average of hourly precipitation in the month of January from
1928 to 1997 recorded in 621 different observation stations located in the
contiguous United States. The coordinates of these locations are normalised to
be inside the $2$-dimensional unit square. We first conduct the parameter
estimation using the MLE and the ML$q$E with several different values of the
hyperparameter $q$, and then try to remove the years in the data with the most
outliers and redo the parameter estimation using the MLE, to see whether the
MLE after removing the outliers is closer to the ML$q$E without removing the
outliers for some $q<1$. For the average January data of each year, we first
estimate the empirical variogram of the data, which is an empirical function
that describes the spatial dependence of the random field, using the R package
geoR (ribeiro2007geor); next, we find out the January data of which years
behave differently from other years via a functional boxplot of the
corresponding empirical variogram, using the R package fbplot
(sun2011functional). The resulting functional boxplot for the January data of
all the 70 years is shown in the first sub-figure of Figure 1, where the
dashed-curves represent the outliers identified among the 70 empirical
variograms, and here we identify 5 outliers in total. Additionally, in the
same plot, the black curve shows the median of the data, and the purple area
represents the central 50% region. In the last three sub-figures of Figure 1,
the blue horizontal lines represent the MLE results with the original data,
and the red dashed lines are the MLE results after removing the 5 outliers
from the dataset, while the black curves are the ML$q$E results without
removing the outliers. The blue dots in those sub-figures are the ML$q$E
values for $q=1$, which coincide with the MLE values. We can see that the red
dashed line intersects with the black curve at $q<1$, which means that the
existence of outliers significantly impacts the parameter estimation results,
and the ML$q$E is able to reduce their effect with $q<1$. It suggests that
applying the ML$q$E with some $q<1$ is able to provide parameter estimation
results diminishing the effect of the outliers, with no need to identify and
remove the outliers in the dataset manually beforehand.

In this paper, we first investigate the asymptotic properties of the ML$q$E in
the scenario of Gaussian random field with Matérn kernel. Then, we carry out
simulation studies to investigate the numerical performance of the ML$q$E with
synthetic spatial data, which demonstrate that ML$q$E is able to provide a
more robust and stable parameter estimation in the presence of outliers. This
observation is also corroborated by experiments on real data. Next, we develop
methodology to choose the optimal value of the hyper-parameter $q$ for the
ML$q$E using its asymptotic variance. Since spatial data are often large-scale
in practice, we implement our method using the software ExaGeoStat, which
provides a high performance computing framework for spatial and spatio-
temporal data, making our implementation capable of dealing with large-scale
data efficiently.

The remainder of this article is organised as follows. Section 2 describes the
proposed methodology and provides some of its theoretical properties on
Gaussian random fields, with some of the corresponding proofs placed in
Appendix A. We also show here our procedure of tuning the hyper-parameter $q$
for the ML$q$E, as well as our algorithm of computing the l$q$-likelihood. In
Section 3, we present the data generation method and numerical results from
simulation studies for synthetic data, and compare the behaviour of the ML$q$E
and the MLE on both clean data and contaminated data, while some additional
experimental results can be found in Appendix B. The experiments on our
mechanism of choosing the optimal value of $q$ is included in Section 3 as
well. Section 4 contains a real application of the ML$q$E on precipitation
data of the United States. The conclusion and discussion are presented in
Section 5.

##  2 Method

###  2.1 Maximum L$q$-Likelihood Estimator

Let $f_{1}$ and $f_{2}$ be two density functions, and $X$ be a random variable
whose density exists. The $q$-entropy (tsallis1988possible) of $f_{2}$ with
respect to $f_{1}$ is defined as

$\mathcal{H}_{q}(f_{1},f_{2})=-\mathbb{E}_{f_{1}}\\{L_{q}[f_{2}(X)]\\},\quad
q>0,$ (2.1)

where

$L_{q}(u)=\begin{cases}\log u,\quad\text{if}\quad q=1;\\\
\left(u^{1-q}-1\right)/(1-q),\quad\text{otherwise}.\end{cases}$ (2.2)

The Maximum L$q$-Likelihood Estimator (ML$q$E) is defined through the
empirical $q$-entropy, where we replace $f_{1}$ and $f_{2}$ in (2.1) by the
empirical density and the density of the parametric model respectively.
Consider an i.i.d. sample $X_{1},\dots,X_{n}$ be from the distribution
$g(x;\bm{\theta}_{0})$ with parameter
$\bm{\theta}_{0}\in\bm{\Theta}\in\mathbb{R}^{p}$ for some integer $p\geq 1$,
the ML$q$E for the true parameter $\bm{\theta}_{0}$ is formulated as

$\tilde{\bm{\theta}}=\operatorname*{arg\,max}_{\bm{\theta}\in\bm{\Theta}}\sum_{i=1}^{n}L_{q}\\{g(X_{i};\bm{\theta})\\},\quad
0<q\leq 1.$ (2.3)

From this definition we can see that, when $q=1$, the ML$q$E is exactly the
same as the MLE, as in this case, the function L$q$ is nothing but the log
function. In addition, it can be shown that, when $q\to 1$, if the ML$q$E
$\tilde{\bm{\theta}}$ exists, it approaches the MLE $\hat{\bm{\theta}}$.

###  2.2 ML$q$E for Gaussian Random Fields

The main target of this work is to apply the ML$q$E on spatial data with
replicates. We consider $m$ sets of independent realizations of an identical
zero-mean Gaussian random field, with $m>1$, denoted by
$\\{\mathbf{Z}_{i}(\mathbf{s)}:\mathbf{s}\in\mathbb{R}^{d}\\}$, where
$i=1,\ldots,m$ and $d\in\mathbb{Z}^{+}$. Here $d$ denotes the dimension of the
locations. For most geographical and environmental datasets, we have $d=2$,
and this will be the setting of all the datasets that we use in the numerical
experiments of this work as well. The spatial covariance function of each
$\mathbf{Z}_{i}$ is given by the Matérn covariance function shown in (1.1). We
consider the case that all of the $m$ realizations of the random fields,
denoted by $\\{\mathbf{Z}_{i},i=1,\dots,m\\}$, are recorded on a fixed set of
$n$ locations, where $n>1$. Therefore, each $\mathbf{Z}_{i}$ is an
$n$-dimensional vector with identical zero-mean multivariate Gaussian
distribution with common density

$\displaystyle
f(\mathbf{z};\bm{\theta})=\frac{1}{\left|2\pi\mathbf{\Sigma}_{\mathcal{M}}\right|^{1/2}}\exp\left(-\frac{1}{2}\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{z}\right),$
(2.4)

where $\mathbf{\Sigma}_{\mathcal{M}}$ is the $n\times n$ covariance matrix
obtained from the Matern covariance function $\mathcal{M}(h,\bm{\theta})$ in
(1.1) based on the $n$ locations.

It follows that the ML$q$E $\tilde{\bm{\theta}}$ (2.3) for
$\\{\mathbf{Z}_{i},i=1,\dots,m\\}$ can be expressed as

$\tilde{\bm{\theta}}=\operatorname*{arg\,max}_{\bm{\theta}\in\bm{\Theta}}\sum_{i=1}^{m}L_{q}\\{f(\mathbf{Z}_{i};\bm{\theta})\\},\quad
0<q\leq 1,$

where the function $\kappa$ is as in (2.4) and L$q$ is as in (2.2).

Define

$\displaystyle\mathbf{U}\left(\mathbf{z};\bm{\theta}\right)=\frac{\partial}{\partial\bm{\theta}}f\left(\mathbf{z};\bm{\theta}\right),$

then for $\\{\mathbf{Z}_{i},i=1,\dots,m\\}$, the ML$q$E $\tilde{\bm{\theta}}$
can alternatively be represented by the solution of the following equation:

$\sum_{i=1}^{m}\mathbf{U}(\mathbf{Z}_{i};\tilde{\bm{\theta}})f(\mathbf{Z}_{i};\tilde{\bm{\theta}})^{1-q}=\mathbf{0},\quad
0<q\leq 1.$ (2.5)

###  2.3 Asymptotic Properties

In this subsection, we aim to prove the consistency and derive the asymptotic
variance of the ML$q$E for Gaussian random fields with Matérn covariance.

Note that the density function (2.4) can be rewritten as

$\displaystyle
f(\mathbf{z};\bm{\theta})=\exp\left[\left\\{\bm{\eta}\left(\bm{\theta}\right)\right\\}^{\top}\mathbf{b}\left(\mathbf{z}\right)-A\left(\bm{\theta}\right)\right],$

where we define

$\displaystyle\mathbf{b}\left(\mathbf{z}\right)=\begin{bmatrix}\mathbf{z}\\\
\text{vec}\left({\mathbf{z}\mathbf{z}^{\top}}\right)\end{bmatrix},\qquad\bm{\eta}\left(\bm{\theta}\right)=\begin{bmatrix}\mathbf{0}_{n}\\\
\text{vec}\left(-\frac{1}{2}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\end{bmatrix},$

and $A\left(\bm{\theta}\right)$ is a normalizing constant. Recall that
according to (1.1), the dimension of $\bm{\theta}$ is 3, which is smaller than
the dimension of the vector $\bm{\eta}\left(\bm{\theta}\right)$ for any
integer $n>1$. Therefore, $f(\mathbf{z};\bm{\theta})$ belongs to a curved
exponential family, see p. 25 and Note 10.6 in p. 79 in lehmann2006theory.

In section 3 of lq, the asymptotic distribution theory for MLqE was derived
for full-rank exponential families. However, the parameters for curved
exponential families estimated via the MLqE also follow similar asymptotic
normality as derived in Theorem 3.2 in lq, and further generalised in Theorem
4.2 in the same article. Making use of these results, we have the following
lemma and theorem.

######  Lemma 2.1.

Let $q_{m}>0$ satisfy that $q_{m}\to 1$ as $m\to\infty$, and the underlying
parameter value $\bm{\theta}_{0}$ of the sample is an interior point of the
parameter space $\bm{\Theta}$, which is compact. Then, the probability that
the Equation (2.5) has a unique solution $\tilde{\bm{\theta}}$ converges to 1
as $n\to\infty$, and $\tilde{\bm{\theta}}\stackrel{{\scriptstyle
P}}{{\longrightarrow}}\bm{\theta}_{0}$ as $n\to\infty$.

###### Proof.

This lemma follows directly from Theorem 3.2 in lq. ∎

Next, the following theorem states asymptotic Gaussianity of the MLqE along
with the form of the asymptotic dispersion matrix. Before stating the theorem,
we first define some necessary quantities for it.

Define

$\displaystyle\mathbf{U}^{*}_{n}(\mathbf{Z},\bm{\theta},q)=$
$\displaystyle\left(\frac{1}{(2\pi)^{\frac{n}{2}}|\mathbf{\Sigma}_{\mathcal{M}}|^{\frac{1}{2}}}\right)\exp\left(-\frac{1-q}{2}\left(\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}\right)\right)$
$\displaystyle\times\left[\frac{1}{2}\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}-\frac{1}{2}\text{tr}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)\right],$
$\displaystyle\mathbf{V}^{*}_{n}(\mathbf{Z},\bm{\theta},q)=$
$\displaystyle(1-q)\left(\frac{1}{(2\pi)^{\frac{n}{2}}|\mathbf{\Sigma}_{\mathcal{M}}|^{\frac{1}{2}}}\right)\exp\left(-\frac{1-q}{2}\left(\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}\right)\right)$
$\displaystyle\times\left[\frac{1}{2}\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}-\frac{1}{2}\text{tr}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)\right]^{2}$
$\displaystyle+\left(\frac{1}{(2\pi)^{\frac{n}{2}}|\mathbf{\Sigma}_{\mathcal{M}}|^{\frac{1}{2}}}\right)\exp\left(-\frac{1-q}{2}\left(\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}\right)\right)$
$\displaystyle\times\left[\frac{1}{2}\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\left[\frac{\partial^{2}\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}^{2}}-2\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\right]\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}\right].$

where $\mathbf{\Sigma}_{\mathcal{M}}$,
$\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}$ and
$\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}$
are functions of $\bm{\theta}$, but we drop the argument $\bm{\theta}$ for
brevity. Define $\bm{\theta}_{m}^{*}$ such that it satisfies

$\displaystyle
E_{\bm{\theta}_{0}}\left[\mathbf{U}^{*}_{n}(\mathbf{Z},\bm{\theta}_{m}^{*},q_{m})\right]=\mathbf{0}.$

Also define

$\displaystyle K_{m}$
$\displaystyle=E_{\bm{\theta}_{0}}\left[\mathbf{U}^{*}_{n}(\mathbf{Z},\bm{\theta}_{m}^{*},q_{m})\mathbf{U}^{*}_{n}(\mathbf{Z},\bm{\theta}_{m}^{*},q_{m})^{\top}\right],$
$\displaystyle J_{m}$
$\displaystyle=E_{\bm{\theta}_{0}}\left[\mathbf{V}^{*}_{n}(\mathbf{Z},\bm{\theta}_{m}^{*},q_{m})\right].$

We have the following theorem:

######  Theorem 2.2.

Let $q_{m}\to 1$ as $m\to\infty$, and the parameter value $\bm{\theta}_{0}$ of
the sample is an interior point of the compact parameter space $\bm{\Theta}$.
Then,

$\displaystyle\sqrt{m}J_{m}^{1/2}K_{m}^{-1/2}J_{m}^{1/2}\left(\tilde{\bm{\theta}}-\bm{\theta}_{m}^{*}\right)\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N_{p}(\mathbf{0}_{n},I_{n})\text{ as }m\to\infty.$

If $\sqrt{m}(q_{m}-1)\to 0$ as $m\to\infty$, we have

$\displaystyle\sqrt{m}J_{m}^{1/2}K_{m}^{-1/2}J_{m}^{1/2}\left(\tilde{\bm{\theta}}-\bm{\theta}_{0}\right)\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N_{p}(\mathbf{0}_{n},I_{n})\text{ as }m\to\infty.$

The proof of this theorem, as well as the expressions of the derivatives
$\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}$ and
$\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}$,
are deferred to Appendix A.

###  2.4 Choice of $q$

Although smaller values of the hyper-parameter $q$ can reduce the influence of
outliers in the data to a larger extent, it also results in much higher
variance and numerical instability, as we will see from the numerical results
in Sections 3.1 and 3.2. Therefore, a mechanism of finding the optimal value
of $q$ is desired. The method we apply to tune the hyper-parameter $q$ is
inspired by the method introduced in ribeiro, with the idea of performing grid
search in a smart way, so as to find an sub-interval of $(0,1]$, such that all
values of $q$ taken from this sub-interval lead to similar estimation results.
For this purpose, we need to find a metric which jointly evaluates the three
parameters that we estimate.

The first option is to follow the idea in ribeiro, which is to standardise the
parameters with their corresponding asymptotic variance, and then sum them up.
Define an ordered grid of values for $q$, denoted by
$q_{0}=1>q_{1}>q_{2}>\dots>q_{K}>0$, and for each $q_{k}$ with $0\leq k\leq
K$, denote the corresponding $p$-dimensional ML$q$E by

$\hat{\bm{\theta}}_{q_{k}}=\left(\hat{\theta}_{q_{k}}^{1},\dots,\hat{\theta}_{q_{k}}^{p}\right)^{\top},$

and define the corresponding vector of standardised estimates
$\mathbf{z}_{q_{k}}$ as

$\mathbf{z}_{q_{k}}=\left(\frac{\hat{\theta}_{q_{k}}^{1}}{\sqrt{m}\,\text{se}\left(\hat{\theta}_{q_{k}}^{1}\right)},\dots,\frac{\hat{\theta}_{q_{k}}^{p}}{\sqrt{m}\,\text{se}\left(\hat{\theta}_{q_{k}}^{p}\right)}\right)^{\top},$

where $m$ is the number of replicates of the spatial data, and
$\text{se}(\hat{\theta})$ denotes the asymptotic standard deviation of
$\hat{\theta}$ derived in Section 2.3, which in practice we calculate in the
following way:

$\text{se}(\hat{\theta}_{q_{k}}^{r})=J_{m}(\hat{\theta}_{q_{k}}^{r})^{-1/2}K_{m}(\hat{\theta}_{q_{k}}^{r})^{1/2}J_{m}(\hat{\theta}_{q_{k}}^{r})^{-1/2},$

where

$\displaystyle K_{m}(\hat{\theta}_{q_{k}}^{r})$
$\displaystyle=\mathbf{U}^{*}_{n}(\mathbf{Z},\hat{\theta}_{q_{k}}^{r},q_{k})\mathbf{U}^{*}_{n}(\mathbf{Z},\hat{\theta}_{q_{k}}^{r},q_{k})^{\top},$
$\displaystyle J_{m}(\hat{\theta}_{q_{k}}^{r})$
$\displaystyle=\mathbf{V}^{*}_{n}(\mathbf{Z},\hat{\theta}_{q_{k}}^{r},q_{k}),$

for $1\leq r\leq p$.

Furthermore, we define the standardised quadratic variation (SQV) for each
$q_{k}$ with $1\leq k\leq K$ as

$\text{SQV}_{q_{k}}=||\mathbf{z}_{q_{k-1}}-\mathbf{z}_{q_{k}}||/p,$

and choose the value of $q$ based on this value.

The corresponding alrogithm of finding the optimal value of $q$, denoted by
$q^{*}$, via the SQV, is shown in Algorithm 1. This algorithm basically
follows what is proposed in ribeiro, which aims to find an optimal $q^{*}$
such that the ML$q$E is stabilised in a neighborhood of $q^{*}$. However, if
no optimal value of $q^{*}$ can be found in this way, the algorithm still
returns $q^{*}=1$: we may not be able to obtain robustness in this case, so we
still choose the MLE to ensure stability.

Input : An ordered grid for $q$: $q_{0}=1>q_{1}>\dots>q_{K}=q_{min}>0$, a
threshold for the SQV: $L>0$, a threshold for the difference between values of
$q$: $\epsilon>0$

while  _$q_{0}-q_{min} >\epsilon$,_ do

Initialise $q^{*}\leftarrow 1$;

Calculate $\text{SQV}_{q_{k}}$ for each $k=1,\dots,K$;

if  _$\forall k$ , $\text{SQV}_{q_{k}}<L$_ then

$q^{*}\leftarrow q_{0}$;

break;

end if

Let $k^{*}$ be the largest integer among $1,\dots,K$ such that
$\text{SQV}_{q_{k^{*}}}\geq L$;

Define a new equally spaced grid for $q$:
$q_{0}=q_{k^{*}}>q_{1}>\dots>q_{K}=q_{min}>0$;

end while

Output : $q^{*}$

Algorithm 1 Mechanism of tuning the hyper-parameter $q$ using the $SQV$

However, directly summing up the standardised values of the parameters may not
be the best way to evaluate them, since the resulting quantity does not have
any statistical meaning. To solve this issue, another option we consider is to
make use of the function $\kappa$ defined as

$\kappa(\sigma^{2},\beta,\nu)=\sigma^{2}\beta^{-2\nu},$ (2.6)

which is proved in zhang to be able to identify Matérn covariance parameters
under infill asymptotics. Therefore, compared with the previous option using
the SQV, the function $\kappa$ is able to provide an evaluation on the fitness
of the model with the estimated parameters. Moreover, the computation time
using the function $\kappa$ is much faster, as computing the SQV requires
computing the asymptotic standard errors of the parameters, which is extremely
time consuming since it involves a lot of complicated large matrix operations.

Similar to the previous option, for the parameter estimation results using
each $q_{k}$ with $0\leq k\leq K$, we calculate the corresponding value of
function $\kappa$, denoted by $\kappa_{q_{k}}$. The algorithm of finding the
optimal value of $q$, denoted by $q^{*}$, via the function $\kappa$, is shown
in Algorithm 2.

Input : An ordered grid for $q$: $q_{0}=1>q_{1}>\dots>q_{K}=q_{min}>0$, a
threshold coefficient: $L>0$, a threshold for the difference between values of
$q$: $\epsilon>0$

while  _$q_{0}-q_{min} >\epsilon$,_ do

Initialise $q^{*}\leftarrow 1$;

Calculate $\kappa_{q_{k}}$ for each $k=0,\dots,K$;

Calculate $d\kappa_{q_{k}}=|\kappa_{q_{k-1}}/\kappa_{q_{k}}-1|$ for each
$k=1,\dots,K$;

if  _$\max_{k}\\{d\kappa_{q_{k}}\\} <L*\min_{k}\\{d\kappa_{q_{k}}\\}$_ then

$q^{*}\leftarrow q_{0}$;

break;

end if

Let $k^{*}$ be the largest integer among $1,\dots,K$ such that
$d\kappa_{q_{k^{*}}}\geq L*\min_{k}\\{d\kappa_{q_{k}}\\}$;

Define a new equally spaced grid for $q$:
$q_{0}=q_{k^{*}}>q_{1}>\dots>q_{K}=q_{min}>0$;

end while

Output : $q^{*}$

Algorithm 2 Mechanism of tuning the hyper-parameter $q$ using the function
$\kappa$

Apart from the evaluating metric, another change that we make here is that,
instead of the fixed threshold in Algorithm 1, we apply a floating threshold
based on the values of $\kappa_{q_{k}}$, with a threshold coefficient $L$. In
our simulation studies in Section 3, we set the value of the coefficient $L$
equal to $4$.

In general, the results given by Algorithms 1 and 2 are quite similar, while
the results from Algorithm 2 are more stable than from Algorithm 1. However,
the computation using Algorithm 2 is much faster than Algorithm 1 due to its
simplicity. Therefore, in Sections 3 and 4, in which we conduct numerical
experiments on both simulated and real world datasets, we tune the hyper-
parameter $q$ using Algorithm 2.

###  2.5 Computation

The computation of the ML$q$E is implemented in the software ExaGeoStat
(exageostat). To make the computation easier and more compatible with the
environment in ExaGeoStat, while evaluating the L$q$-likelihood, we evaluate
the log-likelihood first, and perform an additional transformation on the log-
likelihood if $q\neq 1$.

From (2.4), the log-likelihood for a realization $\mathbf{Z}$ of length $n$
from a Gaussian random field with covariance matrix
$\mathbf{\Sigma}_{\mathcal{M}}$ parametrised by $\bm{\theta}$ can be expressed
as

$l(\mathbf{Z};\bm{\theta})=-\frac{n}{2}\log(2\pi)-\frac{1}{2}\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}-\frac{1}{2}\log|\mathbf{\Sigma}_{\mathcal{M}}|.$

Furthermore, from (2.2), the expression of the L$q$-likelihood in this case is

$L_{q}(\mathbf{Z};\bm{\theta})=\begin{cases}l(\mathbf{Z};\bm{\theta}),\quad\text{if}\,\,q=1;\\\
\left\\{\exp[l(\mathbf{Z};\bm{\theta})\times(1-q)]-1\right\\}/(1-q),\quad\text{otherwise}.\end{cases}$

Here in Algorithm 3, we describe how we evaluate the L$q$-likelihood for an
observation vector $\mathbf{Z}$, from a Gaussian random field with Matérn
covariance parametrised by $\bm{\theta}$.

Input : Observation vector $\mathbf{Z}$ of length $n$, parameter $\bm{\theta}$
(initial value or the value from the previous iteration of optimization), the
$n$ locations

Calculate the Matérn covariance matrix $\bm{\Sigma}$ using $\bm{\theta}$ and
the locations based on (1.1);

$\mathbf{L}\mathbf{L}^{\top}=\bm{\Sigma}$: Cholesky factorization;

Calculate $\mathbf{Z}_{new}$ by solving the linear system
$\mathbf{L}\mathbf{L}^{\top}\mathbf{Z}_{new}=\mathbf{Z}$;

$logdet\leftarrow\log|\bm{\Sigma}|$;

$prod\leftarrow\mathbf{Z}_{new}^{\top}\mathbf{Z}_{new}$;

$lq\leftarrow-0.5\cdot prod-0.5\cdot logdet-0.5\cdot n\cdot\log(2\pi)$;

if  _$q\neq 1$_ then

$lq\leftarrow\\{\exp[lq\times(1-q)]-1\\}/(1-q)$;

end if

Output : The L$q$ likelihood value $lq$

Algorithm 3 Calculation of the L$q$-likelihood

In practice, since $q$ is fixed throughout one single optimization process,
when $q\neq 1$, to simplify the calculation in the final step of the
algorithm, we only need to let $lq\leftarrow\exp[lq\times(1-q)]$ instead of
evaluating the full expression. Moreover, since the value of the log-
likelihood can be very small, the term $\exp[lq\times(1-q)]$ is likely to be
exactly $0$ on most computing systems when the number of locations $n$ gets
large. To solve this issue, again in the final step of the algorithm, when
$q\neq 1$, we let $lq\leftarrow\exp[(lq+n)\times(1-q)]$, which means we
multiply the the L$q$-likelihood by the constant $\exp[n\times(1-q)]$, so as
to make the computation possible.

##  3 Simulation Study

In this section, we conduct simulation studies on the proposed method, in
which we use the software ExaGeoStat to both generate the data and estimate
the parameters. The optimization algorithm we use for the parameter estimation
is BOBYQA (powell2009bobyqa), which is a bound constrained algorithm without
using derivatives embedded in ExaGeoStat.

In Section 2.5, we already showed how we evaluate the L$q$-likelihood for the
optimization. In Algorithm 4, we describe how we generate synthetic data from
a Gaussian random field with Matérn covariance matrix $\bm{\Sigma}$
parametrised by $\bm{\theta}$, with $n$ locations and $m$ replicates, and
level of contamination $r$. If $r=0$, then we are generating clean data
without outliers, which is used in Section 3.1; if $r>0$, which means that the
synthetic data needs to be contaminated by outliers and is used in Section
3.2, then we generate the clean data first, and add noises to part of the data
afterwards.

Input : Fixed locations, true parameter $\bm{\theta}$, number of locations
$n$, number of replicates $m$, level of contamination $r$ with $0\leq r<1$

Calculate the Matérn covariance matrix $\bm{\Sigma}$ using $\bm{\theta}$ and
the locations, based on (1.1);

$\mathbf{L}\mathbf{L}^{\top}=\bm{\Sigma}$: Cholesky factorization;

for  _i in $1:m$_ do

Generate vector $\mathbf{e}_{i}$ of length $n$ from i.i.d. standard normal
distribution;

$\mathbf{Z}_{i}\leftarrow\mathbf{L}\cdot\mathbf{e}_{i}$;

if  _$r >0$_ then

Generate $r_{i}$ from continuous uniform distribution between $0$ and $1$;

if  _$r_{i} <r$_ then

Generate random noise vector $\mathbf{z}_{i}$ of length $n$ from i.i.d. normal
(or other reasonable distributions);

$\mathbf{Z}_{i}\leftarrow\mathbf{Z}_{i}+\mathbf{z}_{i}$;

end if

end if

end for

Output : Realizations $\\{\mathbf{Z}_{i},i=1,\dots,m\\}$

Algorithm 4 Generating synthetic data for simulation experiments

In each of the experiments, we show the estimation results of the variance,
range and smoothness parameters in each setting. In addition, we plot the
estimated value of the function $\kappa$ in (2.6), to evaluate the fitness of
the model with the estimated parameters. Throughout this section and Appendix
B, in the figures where we show our experimental results, the “function
$\kappa$” are referring to the function $\kappa$. Moreover, for all the
simulation experiments in this section, we tune the value of $q$ using
Algorithm 2.

###  3.1 Clean Data

First of all, we conduct some experiments to compare the performance of ML$q$E
with MLE when there are no outliers in the simulated data.

In Figure 2, we show the experimental results with simulated data from Matérn
covariance function with $\sigma^{2}=1,\beta=0.1,\nu=0.5$, which is
essentially the exponential covariance function. Here we simulate the data
with $n=1,600$ locations and $m=100$ replicates, and the experiment is
repeated $100$ times on $100$ different datasets generated from the same
random field to make the boxplots. We use red horizontal lines to indicate the
the true values of the parameters and the function $\kappa$ in (2.6), and blue
vertical lines correspond to the value of $q$ leading to the smallest mean
squared error (MSE). We can see that, in this case, the MLE is giving the best
performance in terms of both consistency and efficiency. We also notice that
when $q$ is close to $1$, the ML$q$E is converging to MLE as expected.

However, in real applications, we do not have the luxury to conduct parameter
estimation on different datasets from the same random field in order to find
the value of $q$ with the smallest MSE. Therefore, we also conduct choice-
of-$q$ experiments using Algorithm 2 for each of the $100$ datasets
individually, to examine whether the optimal value of $q$ can also be found in
this way. In the last column of the boxplot, we present the parameter
estimation results using the values of $q$ chosen by Algorithm 2, which have
slightly larger bias than the MLE results with $q=1$. This is also what we
expect as the compromise of using a robust method, that it results in larger
variance than the MLE when applied on clean data.

Figure 2: The MLE and ML$q$E estimation results with no outlier in the data.
The red horizontal lines correspond to the true values of the parameters or
the function $\kappa$ in (2.6), and the blue vertical lines correspond to the
value of $q$ leading to the smallest mean squared error (MSE). The last column
with the label “s” indicates the estimation results using the $q$ values
selected using Algorithm 2.

###  3.2 Contaminated Data

Next, we test the performance of the ML$q$E and compare it with the MLE when
there are different types of outliers in the simulated data. For all the
results shown here in Figures 4 to 6, the data are generated from Gaussian
random field with Matérn covariance function with
$\sigma^{2}=1,\beta=0.1,\nu=0.5$. Still, we simulate the data with $n=1,600$
locations and $m=100$ replicates, and the experiment is repeated $100$ times
on $100$ different datasets generated from the same random field to make the
boxplots. Experimental results with some other settings can be found in
Appendix B. The mechanism of contaminating the data is as described in
Algorithm 4.

Figure 3: The MLE and ML$q$E estimation results with the data contaminated by
noises generated from $N(0,1)$ with probability $10\%$. The last column with
the label “s” indicates the estimation results using the $q$ values selected
using Algorithm 2. The values of $q$ we use here are $0.5$, $0.6$, $0.7$,
$0.8$, $0.85$, $0.9$, $0.95$, $0.975$, $0.99$, $0.999$, $0.9999$, $1$. The red
horizontal lines correspond to the true values of the parameters or the
function $\kappa$ in (2.6), and the blue vertical lines correspond to the
value of $q$ leading to the smallest MSE.

Figure 4: The MLE and ML$q$E estimation results with the data contaminated by
noises generated from $N(0,1)$ with probability $20\%$. The last column with
the label “s” indicates the estimation results using the $q$ values selected
using Algorithm 2. The values of $q$ we use here are $0.5$, $0.6$, $0.7$,
$0.8$, $0.85$, $0.9$, $0.95$, $0.975$, $0.99$, $0.999$, $0.9999$, $1$. The red
horizontal lines correspond to the true values of the parameters or the
function $\kappa$ in (2.6), and the blue vertical lines correspond to the
value of $q$ leading to the smallest MSE.

From the results we can see that, with the existence of outliers in the data,
the MLE results deviate significantly from the true values, which is in line
with our expectation. However, as we can see from Figures 4 and 4, when the
outliers are generated by adding noises from $N(0,1)$, even if $q$ is only a
little bit smaller than $q$, then the ML$q$E results are much closer to the
true values. For data with outliers generated by adding noises from $N(0,9)$,
as we show in Figure 6, since the magnitude of outliers are much larger,
smaller values of $q$ are needed to provide better ML$q$E results. In
addition, as before, we use red horizontal lines to indicate the the true
values of the parameters and the function $\kappa$ in (2.6), and blue vertical
lines correspond to the value of $q$ leading to the smallest MSE.

Furthermore, similar to Figure 2, in Figures 4 to 6, we also add one column at
the right of each boxplot, to present the parameter estimation results using
the values of $q$ chosen by Algorithm 2. The choice-of-$q$ experiments are
conducted on each of the $100$ datasets individually. We notice that in all of
these three cases, for most of the simulated datasets, the estimation results
using the chosen $q$ are similar to the results using the value of $q$ with
the smallest values of MSE indicated by the blue vertical lines. In Figure 6,
since the outliers in the datasets are generated in a relatively extreme way,
it is more often that the optimal value of $q$ is not chosen, compared with
Figures 4 and 4. Nevertheless, despite those several extreme cases, Algorithm
2 is working well and the resulting parameter estimation results are very
close to the true values.

Figure 5: The MLE and ML$q$E estimation results with the data contaminated by
noises generated from $N(0,9)$ with probability $10\%$. The last column with
the label “s” indicates the estimation results using the $q$ values selected
using Algorithm 2. The values of $q$ we use here are $0.7$, $0.8$, $0.85$,
$0.9$, $0.95$, $0.96$, $0.98$, $0.99$, $0.999$, $0.9999$, $1$. The red
horizontal lines correspond to the true values of the parameters or the
function $\kappa$ in (2.6), and the blue vertical lines correspond to the
value of $q$ leading to the smallest MSE.

Figure 6: The histograms of the values of $q$ chosen by Algorithm 2, with the
same synthetic data as in the experiments shown in Figures 2 to 6. The
$y$-axes indicate the frequency (in percentage) of the selected $q$.

###  3.3 Chosen Values of $q$

In addition to the boxplots, we also present histograms of the chosen values
of $q$ using Algorithm 2 in Figure 6, with the same datasets that we use to
plot Figures 2 to 6.

For the clean data, as expected, $q=1$ is chosen for most of the simulated
datasets, and among the rest, $q>0.99$ is chosen for most cases, with $q<0.99$
chosen for only a very small number of times. For the contaminated data, in
each of the three cases, the most frequently chosen value of $q$ is equal to
the value of $q$ with the smallest MSE, which are $0.99$ when the outliers are
generated by adding noises from $N(0,1)$, and $0.98$ when the outliers are
generated by adding noises from $N(0,9)$. Moreover, recall that the algorithm
to tune $q$ is designed in such a way that it chooses $1$ if the ML$q$E
results are too unstable, which explains why in the most extreme case where
the outliers are generated by adding noises from $N(0,9)$, $q=1$ is chosen for
more datasets than the other cases of contaminated data. Those datasets for
which $q=1$ is chosen also correspond to most of the outlying points in the
columns with the label “s” in Figure 6.

##  4 Precipitation Data Analysis

For the real data analysis, we consider the US precipitation data the same as
which we use for the motivating example in Section 1. To eliminate the
influence of missing data on our experiment, we selected the data from 1928 to
1997 on $n=621$ observation stations whose data do not contain any missing
value for those $m=70$ years. In addition, we extract the data of each month
in the dataset and study them separately. For example, for the experiment
labelled as “Jan”, we use only the precipitation data of January from 1928 to
1997 in all the observation stations. In this way, it makes more sense to
assume that all the replicates (i.e., data from different years) in each
experiment are from the same distribution. Moreover, to remove the trend in
the data, the data for each year are centralised separately.

In Figure 7, we show the experimental results using the real data. Within each
of these 3 figures, in the first 3 sub-figures, we show the estimation results
using ML$q$E with different values of $q$, while it degenerates to the MLE
when $q=1$; in the last sub-figure, we show the estimated value of the
function shown in (2.6). The values of $q$ that we use here are $0.9$,
$0.925$, $0.95$, $0.96$, $0.97$, $0.98$, $0.99$, $0.999$, $0.9999$, and $1$.
The dot on each curve shows the optimal value of $q$ selected for the
corresponding data using Algorithm 2.

Figure 7: The MLE and ML$q$E estimation results with different values of $q$,
for the precipitation data of different months. Each row shows the plot of the
estimation of the same parameter or the function $\kappa$ (2.6), and each
column shows the results for the same months. The legends in the second row of
each column apply to all sub-figures in the same column. The dots on the
curves indicate the values of $q$ selected with Algorithm 2 and their
corresponding estimation results. The abbreviations in the titles, “DJF”,
“MAM”, “JJA” and “SON”, represent the months in each season, namely December,
January and February for Winter, March, April and May for Spring, June, July
and August for Summer, and September, October and November for Autumn

From the results shown in Figure 7, we can see that for all the 12 months, the
ML$q$E results with most values of $q<1$ are different from the corresponding
MLE results. For the data of October, the selected optimal value of $q$ is
$1$, because the ML$q$E results, especially for the range parameter, are not
able to stabilise as $q$ decreases, and the MLE is considered as the optimal
choice to ensure that we still have a “not too bad” estimation result. For the
data of each of the other months, The optimal values of $q$ are all selected
to be a point such that the curve of the estimation results becomes relatively
flat in a neighborhood of it, which means that stable ML$q$E results can be
achieved using the $q$ values closed to this value.

##  5 Conclusion

To conclude, both simulation and real data experiments have shown that, the
ML$q$E can indeed outperform the MLE with the existence of outliers in the
data. With a suitable choice of $q$, the parameter estimation results via the
ML$q$E exhibit extraordinary low bias and high efficiency, compared with the
MLE, when the data is contaminated. On the theory side, the ML$q$E holds
desirable asymptotic properties as well. The mechanism of tuning the hyper-
parameter $q$ for the ML$q$E that we develop is also proved to be working very
well by our experiments. Moreover, our method is able to efficiently deal with
data with much larger scale than in literature, with the help of the high
performance computing framework ExaGeoStat, which is extremely important for
spatial data analysis.

A limitation of the current work is that the current method is only able to
provide robust estimation for spatial data with multiple replicates, since it
can only treat some of the replicates as outliers as a whole. One interesting
follow up work would be investigating on how to adjust this method so that it
can be applied on spatial data with only one replicate as well. Another
possible extension of this work would be integrating the proposed method with
high performance approximation methods for numerical linear algebra, such as
the tile low-rank algorithm and the mixed-precision method, so that the ML$q$E
can be further applied on even larger-scale datasets efficiently.

##  Appendix A Proofs

###  A.1 Proof of Theorem 2.2

###### Proof.

Note that

$\displaystyle\log f(\mathbf{z};\bm{\theta})$
$\displaystyle=-\frac{n}{2}\log(2\pi)-\frac{1}{2}\log|\mathbf{\Sigma}_{\mathcal{M}}|-\frac{1}{2}\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{z},$
$\displaystyle\frac{\partial}{\partial\bm{\theta}}\log|\mathbf{\Sigma}_{\mathcal{M}}|$
$\displaystyle=\text{tr}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right),$
$\displaystyle\frac{\partial}{\partial\bm{\theta}}\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}\mathbf{z}$
$\displaystyle=\mathbf{z}^{\top}\left(\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\mathbf{z}=-\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{z}.$

Hence, we have

$\displaystyle\frac{\partial}{\partial\bm{\theta}}\log
f(\mathbf{z};\bm{\theta})$
$\displaystyle=-\frac{1}{2}\text{tr}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)+\frac{1}{2}\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{z}.$

Therefore,

$\displaystyle\frac{\partial}{\partial\bm{\theta}}L_{q}[f(\mathbf{Z};\bm{\theta})]=\frac{\partial}{\partial\bm{\theta}}\frac{f^{1-q}(\mathbf{Z};\bm{\theta})-1}{1-q}=f^{1-q}(\mathbf{Z};\bm{\theta})\frac{\partial}{\partial\bm{\theta}}\log
f(\mathbf{Z};\bm{\theta})$ $\displaystyle=$
$\displaystyle\left(\frac{1}{(2\pi)^{\frac{n}{2}}|\mathbf{\Sigma}_{\mathcal{M}}|^{\frac{1}{2}}}\right)\exp\left(-\frac{1-q}{2}\left(\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}\right)\right)\left[\frac{1}{2}\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}-\frac{1}{2}\text{tr}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)\right].$
(A.1)

From (A.1) we get that the expression of
$\mathbf{U}^{*}_{n}(\mathbf{Z},\bm{\theta},q)$ defined here matches the
quantity $\mathbf{U}^{*}(X,\theta,q)$ defined in (2.6) in lq. Therefore,
$\bm{\theta}_{m}^{*}$ defined here corresponds to $\theta_{n}^{*}$ defined in
(3.2) in lq.

Next,

$\displaystyle\frac{\partial^{2}}{\partial\bm{\theta}^{2}}L_{q}[f(\mathbf{Z};\bm{\theta})]$
$\displaystyle=\left(\frac{\partial}{\partial\bm{\theta}}f^{1-q}(\mathbf{Z};\bm{\theta})\right)\frac{\partial}{\partial\bm{\theta}}\log
f(\mathbf{Z};\bm{\theta})+f^{1-q}(\mathbf{Z};\bm{\theta})\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\log
f(\mathbf{Z};\bm{\theta})$
$\displaystyle=(1-q)f^{1-q}(\mathbf{Z};\bm{\theta})\left(\frac{\partial}{\partial\bm{\theta}}\log
f(\mathbf{Z};\bm{\theta})\right)^{2}+f^{1-q}(\mathbf{Z};\bm{\theta})\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\log
f(\mathbf{Z};\bm{\theta}).$

We proceed to derive the expression of
$\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\log f(\mathbf{Z};\bm{\theta})$.
Note that

$\displaystyle\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\log|\mathbf{\Sigma}_{\mathcal{M}}|$
$\displaystyle=\frac{\partial}{\partial\bm{\theta}}\text{tr}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)=\text{tr}\left(\frac{\partial}{\partial\bm{\theta}}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)\right)$
$\displaystyle=\text{tr}\left(\left(\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\left(\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)+\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}\right),$
$\displaystyle\frac{\partial}{\partial\bm{\theta}}\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}\mathbf{z}$
$\displaystyle=\mathbf{z}^{\top}\left(\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\mathbf{z}=-\mathbf{z}^{\top}\left(\frac{\partial}{\partial\bm{\theta}}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\right)\mathbf{z},$
$\displaystyle\frac{\partial}{\partial\bm{\theta}}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)$
$\displaystyle=\left(\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}+\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\left(\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)$
$\displaystyle=\left(\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}+\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial^{2}\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}+\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\left(\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)$
$\displaystyle=-2\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}+\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial^{2}\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}$
$\displaystyle=\mathbf{\Sigma}_{\mathcal{M}}^{-1}\left[\frac{\partial^{2}\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}^{2}}-2\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\right]\mathbf{\Sigma}_{\mathcal{M}}^{-1},$
$\displaystyle\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}\mathbf{z}$
$\displaystyle=-\mathbf{z}^{\top}\left(\frac{\partial}{\partial\bm{\theta}}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\right)\mathbf{z}$
$\displaystyle=\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\left[-\frac{\partial^{2}\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}^{2}}+2\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\right]\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{z}.$

Hence,

$\displaystyle\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\log
f(\mathbf{z};\bm{\theta})=$
$\displaystyle-\frac{1}{2}\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\log|\mathbf{\Sigma}_{\mathcal{M}}|-\frac{1}{2}\frac{\partial}{\partial\bm{\theta}}\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}\mathbf{z}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\text{tr}\left(\left(\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\right)\left(\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)+\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}\right)$
$\displaystyle-\frac{1}{2}\mathbf{z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\left[-\frac{\partial^{2}\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}^{2}}+2\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\right]\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{z},$

and it follows that

$\displaystyle\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\log
f(\mathbf{Z};\bm{\theta})$ $\displaystyle=$
$\displaystyle(1-q)f^{1-q}(\mathbf{Z};\bm{\theta})\left(\frac{\partial}{\partial\bm{\theta}}\log
f(\mathbf{Z};\bm{\theta})\right)^{2}+f^{1-q}(\mathbf{Z};\bm{\theta})\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\log
f(\mathbf{Z};\bm{\theta})$ $\displaystyle=$
$\displaystyle(1-q)\left(\frac{1}{(2\pi)^{\frac{n}{2}}|\mathbf{\Sigma}_{\mathcal{M}}|^{\frac{1}{2}}}\right)\exp\left(-\frac{1-q}{2}\left(\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}\right)\right)\left[\frac{1}{2}\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}-\frac{1}{2}\text{tr}\left(\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}\right)\right]^{2}$
$\displaystyle+\left(\frac{1}{(2\pi)^{\frac{n}{2}}|\mathbf{\Sigma}_{\mathcal{M}}|^{\frac{1}{2}}}\right)\exp\left(-\frac{1-q}{2}\left(\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}\right)\right)\left[\frac{1}{2}\mathbf{Z}^{\top}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\left[\frac{\partial^{2}\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}^{2}}-2\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}^{-1}\frac{\partial\mathbf{\Sigma}_{\mathcal{M}}}{\partial\bm{\theta}}\right]\mathbf{\Sigma}_{\mathcal{M}}^{-1}\mathbf{Z}\right].$

Therefore, the expressions of $K_{m}$ and $J_{m}$ defined in the statement of
the theorem matches with $K_{n}$ and $J_{n}$ defined in (3.4) and (3.5) in lq,
respectively, which completes the proof with reference to Theorem 3.2 and
Theorem 4.3 in lq.

∎

###  A.2 Expressions of Derivatives

For practical use of Theorem 2.2, we need to evaluate
$\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}$ and
$\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}$.
Note that $\frac{\partial}{\partial\bm{\theta}}\mathbf{\Sigma}_{\mathcal{M}}$
and
$\frac{\partial^{2}}{\partial\bm{\theta}^{2}}\mathbf{\Sigma}_{\mathcal{M}}$
are arrays consisting of the second order derivatives of the terms of
$\mathbf{\Sigma}_{\mathcal{M}}$ with respect to $\bm{\theta}$. The terms of
$\mathbf{\Sigma}_{\mathcal{M}}$ are of the form $\mathcal{M}(h,\theta)$, as
defined below.

Denote

$\mathcal{M}(h,\theta)=\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right).$

Then we have

$\displaystyle\frac{\partial\mathcal{M}(h,\theta)}{\partial\sigma^{2}}=$
$\displaystyle\frac{1}{\Gamma(\nu)2^{\nu-1}}\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right);$
$\displaystyle\frac{\partial^{2}\mathcal{M}(h,\theta)}{\partial(\sigma^{2})^{2}}=$
$\displaystyle 0;$
$\displaystyle\frac{\partial\mathcal{M}(h,\theta)}{\partial\beta}=$
$\displaystyle\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left\\{-\frac{\nu}{\beta}\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)-\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}^{\prime}_{\nu}\left(\frac{h}{\beta}\right)\frac{h}{\beta^{2}}\right\\};$
$\displaystyle\frac{\partial^{2}\mathcal{M}(h,\theta)}{\partial\beta^{2}}=$
$\displaystyle\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left\\{\frac{h}{\beta^{2}}\nu\left(\frac{h}{\beta}\right)^{\nu-1}\left[\frac{\nu}{\beta}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)+\mathcal{K}^{\prime}_{\nu}\left(\frac{h}{\beta}\right)\frac{h}{\beta}\right]\right.$
$\displaystyle\left.-\left(\frac{h}{\beta}\right)^{\nu}\left[-\frac{\nu}{\beta^{2}}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)-\frac{\nu}{\beta}\frac{h}{\beta^{2}}\mathcal{K}^{\prime}_{\nu}\left(\frac{h}{\beta}\right)-\frac{h}{\beta^{2}}\mathcal{K}^{\prime}_{\nu}\left(\frac{h}{\beta}\right)-\frac{h^{2}}{\beta^{3}}\mathcal{K}^{\prime\prime}_{\nu}\left(\frac{h}{\beta}\right)\right]\right\\};$
$\displaystyle\frac{\partial\mathcal{M}(h,\theta)}{\partial\nu}=$
$\displaystyle\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left\\{-\left(\log(2)+\Psi(\nu)\right)\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)+\frac{\partial}{\partial\nu}\left[\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)\right]\right\\};$
$\displaystyle\frac{\partial^{2}\mathcal{M}(h,\theta)}{\partial\nu^{2}}=$
$\displaystyle-\left(\log(2)+\Psi(\nu)\right)\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left\\{-\left(\log(2)+\Psi(\nu)\right)\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)+\frac{\partial}{\partial\nu}\left[\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)\right]\right\\}$
$\displaystyle-\Psi^{\prime}(\nu)\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)$
$\displaystyle+\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left\\{-\left(\log(2)+\Psi(\nu)\right)\frac{\partial}{\partial\nu}\left[\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)\right]+\frac{\partial^{2}}{\partial\nu^{2}}\left[\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)\right]\right\\}.$

Here $\Psi(\cdot)$ represents the digamma function.

For the cross terms:

$\displaystyle\frac{\partial^{2}\mathcal{M}(h,\theta)}{\partial(\sigma^{2})\partial\beta}=$
$\displaystyle\frac{1}{\Gamma(\nu)2^{\nu-1}}\left\\{-\frac{\nu}{\beta}\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)-\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}^{\prime}_{\nu}\left(\frac{h}{\beta}\right)\frac{h}{\beta^{2}}\right\\};$
$\displaystyle\frac{\partial^{2}\mathcal{M}(h,\theta)}{\partial(\sigma^{2})\partial\nu}=$
$\displaystyle\frac{1}{\Gamma(\nu)2^{\nu-1}}\left\\{-\left(\log(2)+\Psi(\nu)\right)\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)+\frac{\partial}{\partial\nu}\left[\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)\right]\right\\};$
$\displaystyle\frac{\partial^{2}\mathcal{M}(h,\theta)}{\partial\beta\partial\nu}=$
$\displaystyle\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left\\{\left(\log(2)+\Psi(\nu)\right)\left[\frac{\nu}{\beta}\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)+\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}^{\prime}_{\nu}\left(\frac{h}{\beta}\right)\frac{h}{\beta^{2}}\right]\right.$
$\displaystyle-\left.\frac{1}{\beta}\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)-\frac{\nu}{\beta}\frac{\partial}{\partial\nu}\left[\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{h}{\beta}\right)\right]-\frac{h}{\beta^{2}}\frac{\partial}{\partial\nu}\left[\left(\frac{h}{\beta}\right)^{\nu}\mathcal{K}^{\prime}_{\nu}\left(\frac{h}{\beta}\right)\right]\right\\}.$

##  Appendix B Supplementary Experimental Results

###  B.1 Choice of $q$ Algorithms

Here we present a comparison between the two methods that we propose in
Section 2.4, Algorithm 1 and 2, using the same data as the experiment shown in
Figure 4. In Figure 8, we show the parameter estimation results with the
values of $q$ chosen by Algorithm 1 (in the columns named “s1”) and Algorithm
2 (in the columns named “s2”), along with the MLE results ($q=1$) and ML$q$E
results with $q=0.99$, which is the $q$ leading to the smallest MSE in Figure
4. Both Algorithms 1 and 2 are able to give quite unbiased estimation results,
while the results in “s2” are more stable than in “s1”, as there are much less
outlying points. Further considering the massive advantage on computation
speed of Algorithm 2 over Algorithm 1, in practice, it would be recommended to
apply Algorithm 2 to tune the hyper-parameter $q$.

Figure 8: The MLE and ML$q$E estimation results with the data contaminated by
noises generated from $N(0,1)$ with probability $10\%$. The last column with
the label “s” indicates the estimation results using the $q$ values selected
using Algorithm 2. The values of $q$ we use here are $0.5$, $0.6$, $0.7$,
$0.8$, $0.85$, $0.9$, $0.95$, $0.975$, $0.99$, $0.999$, $0.9999$, $1$. The red
horizontal lines correspond to the true values of the parameters or the
function $\kappa$ in (2.6), and the blue vertical lines correspond to the
value of $q$ leading to the smallest MSE

###  B.2 Parameter Estimation

In Section 3, we presented the simulation results using spatial data generate
using Matérn kernel with true parameters $\sigma^{2}=1,\beta=0.1,\nu=0.5$. In
this appendix, we present some simulation results with some different values
of $\beta$ and $\nu$, with other settings the same as in Section 3. Generally,
we consider the cases with $\beta=0.03,0.1$, which we refer to as weak and
medium correlation respectively, and $\nu=0.5,1$, which we refer to as rough
and smooth field respectively. For all the experimental results shown here, we
set the number of locations $n=1{,}600$, number of replicates $m=100$, and we
repeat the experiment $100$ times to make the boxplots. All experiments are
conducted using the software ExaGeoStat.

####  B.2.1 Medium-smooth Case

Here we present the results for $\sigma^{2}=1,\beta=0.1,\nu=1$ in Figures 9 to
11.

Figure 9: The MLE and ML$q$E estimation results, with $1,600$ locations and
$100$ replicates, and experiment is repeated $100$ times to make the boxplots.
The replicates in the data are contaminated by noises generated from $N(0,1)$
with probability $10\%$. The true values of the parameters are
$\sigma^{2}=1,\beta=0.1,\nu=1$. The rightmost box represents the MLE, while
the others represent the ML$q$E and the values of the x-axis correspond to the
values of $q$. The values of $q$ that we use here are $0.9$, $0.95$, $0.96$,
$0.97$, $0.98$, $0.99$, $1$. The red horizontal lines correspond to the true
values of the parameters or the function $\kappa$ in (2.6) Figure 10: The MLE
and ML$q$E estimation results, with $1,600$ locations and $100$ replicates,
and experiment is repeated $100$ times to make the boxplots. The replicates in
the data are contaminated by noises generated from $N(0,1)$ with probability
$20\%$. The true values of the parameters are $\sigma^{2}=1,\beta=0.1,\nu=1$.
The rightmost box represents the MLE, while the others represent the ML$q$E
and the values of the x-axis correspond to the values of $q$. The values of
$q$ that we use here are $0.9$, $0.95$, $0.96$, $0.97$, $0.98$, $0.99$, $1$.
The red horizontal lines correspond to the true values of the parameters or
the function $\kappa$ in (2.6) Figure 11: The MLE and ML$q$E estimation
results, with $1,600$ locations and $100$ replicates, and experiment is
repeated $100$ times to make the boxplots. The replicates in the data are
contaminated by noises generated from $N(0,9)$ with probability $10\%$. The
true values of the parameters are $\sigma^{2}=1,\beta=0.1,\nu=1$. The
rightmost box represents the MLE, while the others represent the ML$q$E and
the values of the x-axis correspond to the values of $q$. The values of $q$
that we use here are $0.9$, $0.95$, $0.96$, $0.97$, $0.98$, $0.99$, $1$. The
red horizontal lines correspond to the true values of the parameters or the
function $\kappa$ in (2.6)

####  B.2.2 Weak-rough Case

Here we present the results for $\sigma^{2}=1,\beta=0.03,\nu=0.5$ in Figures
12 to 14.

Figure 12: The MLE and ML$q$E estimation results, with $1,600$ locations and
$100$ replicates, and experiment is repeated $100$ times to make the boxplots.
The replicates in the data are contaminated by noises generated from $N(0,1)$
with probability $10\%$. The true values of the parameters are
$\sigma^{2}=1,\beta=0.03,\nu=0.5$. The rightmost box represents the MLE, while
the others represent the ML$q$E and the values of the x-axis correspond to the
values of $q$. The values of $q$ that we use here are $0.9$, $0.95$, $0.96$,
$0.97$, $0.98$, $0.99$, $1$. The red horizontal lines correspond to the true
values of the parameters or the function $\kappa$ in (2.6) Figure 13: The MLE
and ML$q$E estimation results, with $1,600$ locations and $100$ replicates,
and experiment is repeated $100$ times to make the boxplots. The replicates in
the data are contaminated by noises generated from $N(0,1)$ with probability
$20\%$. The true values of the parameters are
$\sigma^{2}=1,\beta=0.03,\nu=0.5$. The rightmost box represents the MLE, while
the others represent the ML$q$E and the values of the x-axis correspond to the
values of $q$. The values of $q$ that we use here are $0.9$, $0.95$, $0.96$,
$0.97$, $0.98$, $0.99$, $1$. The red horizontal lines correspond to the true
values of the parameters or the function $\kappa$ in (2.6) Figure 14: The MLE
and ML$q$E estimation results, with $1,600$ locations and $100$ replicates,
and experiment is repeated $100$ times to make the boxplots. The replicates in
the data are contaminated by noises generated from $N(0,9)$ with probability
$10\%$. The true values of the parameters are
$\sigma^{2}=1,\beta=0.03,\nu=0.5$. The rightmost box represents the MLE, while
the others represent the ML$q$E and the values of the x-axis correspond to the
values of $q$. The values of $q$ that we use here are $0.9$, $0.95$, $0.96$,
$0.97$, $0.98$, $0.99$, $1$. The red horizontal lines correspond to the true
values of the parameters or the function $\kappa$ in (2.6)

####  B.2.3 Weak-smooth Case

Here we present the results for $\sigma^{2}=1,\beta=0.03,\nu=1$ in Figures 15
and 16.

Figure 15: The MLE and ML$q$E estimation results, with $1,600$ locations and
$100$ replicates, and experiment is repeated $100$ times to make the boxplots.
The replicates in the data are contaminated by noises generated from $N(0,1)$
with probability $10\%$. The true values of the parameters are
$\sigma^{2}=1,\beta=0.03,\nu=1$. The rightmost box represents the MLE, while
the others represent the ML$q$E and the values of the x-axis correspond to the
values of $q$. The values of $q$ that we use here are $0.9$, $0.95$, $0.96$,
$0.97$, $0.98$, $0.99$, $1$. The red horizontal lines correspond to the true
values of the parameters or the function $\kappa$ in (2.6) Figure 16: The MLE
and ML$q$E estimation results, with $1,600$ locations and $100$ replicates,
and experiment is repeated $100$ times to make the boxplots. The replicates in
the data are contaminated by noises generated from $N(0,1)$ with probability
$20\%$. The true values of the parameters are $\sigma^{2}=1,\beta=0.03,\nu=1$.
The rightmost box represents the MLE, while the others represent the ML$q$E
and the values of the x-axis correspond to the values of $q$. The values of
$q$ that we use here are $0.9$, $0.95$, $0.96$, $0.97$, $0.98$, $0.99$, $1$.
The red horizontal lines correspond to the true values of the parameters or
the function $\kappa$ in (2.6)


# Short-Term Meaning Shift: A Distributional Exploration

Marco Del Tredici∗ Raquel Fernández∗ Gemma Boleda†
∗University of Amsterdam †Universitat Pompeu Fabra
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

We present the first exploration of meaning shift over short periods of time
in online communities using distributional representations. We create a small
annotated dataset and use it to assess the performance of a standard model for
meaning shift detection on short-term meaning shift. We find that the model
has problems distinguishing meaning shift from referential phenomena, and
propose a measure of contextual variability to remedy this.

##  1 Introduction

Semantic change has received increasing attention in empirical Computational
Linguistics / NLP in the last few years Tang (2018); Kutuzov et al. (2018).
Almost all studies so far have focused on meaning shift in long periods of
time—decades to centuries. However, the genesis of meaning shift and the
mechanisms that produce it operate at much shorter time spans, ranging from
the online agreement on words’ meaning in dyadic interactions Brennan and
Clark (1996) to the rapid spread of new meanings in relatively small
communities of people in Wenger (1998); Eckert and McConnell-Ginet (1992).
This paper is, to the best of our knowledge, the first exploration of the
latter phenomenon—which we call short-term meaning shift—using distributional
representations.

More concretely, we focus on meaning shift arising within a period of 8 years,
and explore it on data from an online community of speakers, because there the
adoption of new meanings happens at a fast pace Clark (1996); Hasan (2009).
Indeed, short-term shift is usually hard to observe in standard language, such
as the language of books or news, which has been the focus of long-term shift
studies (e.g., Hamilton et al., 2016; Kulkarni et al., 2015), since it takes a
long time for a new meaning to be widely accepted in the standard language.

Our contribution is twofold. First, we create a small dataset of short-term
shift for analysis and evaluation, and qualitatively analyze the types of
meaning shift we find.111Data and code are available at:
https://github.com/marcodel13/Short-term-meaning-shift. This is necessary
because, unlike studies of long-term shift, we cannot rely on material
previously gathered by linguists or lexicographers. Second, we test the
behavior of a standard distributional model of semantic change when applied to
short-term shift. Our results show that this model successfully detects most
shifts in our data, but it overgeneralizes. Specifically, the model gets
confused with contextual changes due to speakers in the community often
talking about particular people and events, which are frequent on short time
spans. We propose to use a measure of contextual variability to remedy this
and showcase its potential to spot false positives of referential nature like
these. We thus make progress in understanding the nature of semantic shift and
towards improving computational models thereof.

##  2 Related Work

Distributional models of semantic change are based on the hypothesis that a
change in context of use mirrors a change in meaning. This in turn stems from
the Distributional Hypothesis, that states that similarity in meaning results
in similarity in context of use Harris (1954). Therefore, all models
(including ours) spot semantic shift as a change in the word representation in
different time periods. Among the most widely used techniques are Latent
Semantic Analysis Sagi et al. (2011); Jatowt and Duh (2014), Topic Modeling
Wijaya and Yeniterzi (2011), classic distributional representations based on
co-occurence matrices of target words and context terms Gulordava and Baroni
(2011). More recently, researchers have used word embeddings computed using
the skip-gram model by Mikolov et al. (2013). Since embeddings computed in
different semantic spaces are not directly comparable, time related
representations are usually made comparable either by aligning different
semantic spaces through a transformation matrix Kulkarni et al. (2015);
Azarbonyad et al. (2017); Hamilton et al. (2016) or by initializing the
embeddings at $t_{i+1}$ using those computed at $t_{i}$ Kim et al. (2014); Del
Tredici et al. (2016); Phillips et al. (2017); Szymanski (2017). We adopt the
latter methodology (see Section 3.2).

Unlike most previous work, we focus on the language of online communities.
Recent studies of this type of language have investigated the spread of new
forms and meanings Del Tredici and Fernández (2017, 2018); Stewart and
Eisenstein (2018), competing lexical variants Rotabi et al. (2017), and the
relation between conventions in a community and the social standing of its
members Danescu-Niculescu-Mizil et al. (2013). None of these works has
analyzed the ability of a distributional model to capture these phenomena,
which is what we do in this paper for short-term meaning shift. Kulkarni et
al. (2015) consider meaning shift in short time periods on Twitter data, but
without providing an analysis of the observed shift nor systematically
assessing the performance of the model, as we do here.

Evaluation of semantic shift is difficult, due to the lack of annotated
datasets Frermann and Lapata (2016). For this reason, even for long-term
shift, evaluation is usually performed by manually inspecting the $n$ words
whose representation changes the most according to the model under
investigation Hamilton et al. (2016); Kim et al. (2014). Our dataset allows
for a more systematic evaluation and analysis, and enables comparison in
future studies.

##  3 Experimental Setup

###  3.1 Data

We exploit user-generated language from an online forum of football fans,
namely, the r/LiverpoolFC subreddit, one of the many communities hosted by the
Reddit platform.222https://www.reddit.com. We downloaded Reddit data using the
Python package Praw: https://pypi.python.org/pypi/praw/. Del Tredici and
Fernández (2018) showed that this subreddit presents many characteristics that
favour the creation and spread of linguistic innovations, such as a topic that
reflects a strong external interest and high density of the connections among
its users. This makes it a good candidate for our investigation. We focus on a
short period of eight years, between 2011 and 2017. In order to enable a
clearer observation of short-term meaning shift, we define two non-consecutive
time bins: the first one ($t_{1}$) contains data from 2011–2013 and the second
one ($t_{2}$) from 2017.333These choices ensure that the samples in these two
time bins are approximately of the same size – see Table 1. The r/LiverpoolFC
subreddit exists since 2009, but very little content was produced in
2009–2010. We also use a large sample of community-independent language for
the initialization of the word vectors, namely, a random crawl from Reddit in
2013. Table 1 shows the size of each sample.

###  3.2 Model

We adopt the model proposed by Kim et al. (2014), a representative method for
computing diachronic meaning shift.444The model was implemented using the
Python package Gensim: https://pypi.python.org/pypi/gensim/. While other
methods might be equally suitable (see Section 2), we expect our results not
to be method-specific, because they concern general properties of short-term
shift, as we show in Sections 4 and 5. In the model by Kim et al. (2014), word
embeddings for the first time bin $t_{1}$ are initialized randomly; then,
given a sequence of time-related samples, embeddings for $t_{i}$ are
initialized using the embeddings of $t_{i-1}$ and further updated. If at
$t_{i}$ the word is used in the same contexts as in $t_{i-1}$, its embedding
will only be marginally updated, whereas a major change in the context of use
will lead to a stronger update of the embedding. The model makes embeddings
across time bins directly comparable.

We implement the following steps: First, we create randomly initialized word
embeddings with the large sample Reddit13 to obtain meaning representations
that are community-independent. We then use these embeddings to initialize
those in LiverpoolFC13, update the vectors on this sample, and thus obtain
embeddings for time $t_{1}$. This step adapts the general embeddings to the
LiverpoolFC community. Finally, we initialize the word embeddings for
LiverpoolFC17 with those of $t_{1}$, train on this sample, and get embeddings
for $t_{2}$.

The vocabulary is defined as the intersection of the vocabularies of the three
samples (Reddit13, LiverpoolFC13, LiverpoolFC17), and includes 157k words. For
Reddit13, we include only words that occur at least 20 times in the sample, so
as to ensure meaningful representations for each word, while for the other two
samples we do not use any frequency threshold: Since the embeddings used for
the initialization of LiverpoolFC13 encode community-independent meanings, if
a word doesn’t occur in LiverpoolFC13 its representation will simply be as in
Reddit13, which reflects the idea that if a word is not used in a community,
then its meaning is not altered within that community. We train with standard
skip-gram parameters Levy et al. (2015): window 5, learning rate 0.01,
embedding dimension 200, hierarchical softmax.

sample | time bin | million tokens
---|---|---
Reddit13 | 2013 |  $\scriptstyle\sim$900
LiverpoolFC13 | 2011–13 |  8.5
LiverpoolFC17 | 2017 | 11.9
Table 1: Time bin and size of the datasets.

###  3.3 Evaluation dataset

Our dataset consists of 97 words from the r/LiverpoolFC subreddit with
annotations by members of the subreddit —that is, community members with
domain knowledge (needed for this task) but no linguistic background.

To ensure that we would get enough cases of semantic shift to enable a
meaningful analysis, we started out from content words that increase their
relative frequency between $t_{1}$ and $t_{2}$.555Frequency increase has been
shown to positively correlate with meaning change Wijaya and Yeniterzi (2011);
Kulkarni et al. (2015); although it is not a necessary condition, it is a
reasonable starting point, as a random selection of words would contain very
few positive examples. Our dataset is thus biased towards precision over
recall. A threshold of 2 standard deviations above the mean yielded $\sim$200
words. The first author manually identified 34 semantic shift candidates among
these words by analyzing their contexts of use in the r/LiverpoolFC data.
Semantic shift is defined here as a change in the ontological type that a word
denotes, which takes place when the word starts to be used to denote an entity
which is different from the one originally denoted and the new use spreads
among the members of a community (see examples in Sec. 4). We added two types
of confounders: 33 words with a significant frequency increase but not marked
as meaning shift candidates, and 33 words with constant frequency between
$t_{1}$ and $t_{2}$, included as a sanity check. All words have absolute
frequency in range [50–500].

The participants were shown the 100 words (in randomized order) together with
randomly chosen contexts of usage from each time period ($\mu$=4.7 contexts
per word) and, for simplicity, were asked to make a binary decision about
whether there was a change in meaning. In order to have the redditors
familiarize with the concept of meaning change, we first provide them with an
intuitive, non-technical definition, and then a set of cases that exemplify
it. The instructions to participants can be found in the project’s GitHub
repository (see footnote 1).

Semantic shift is arguably a graded notion. In line with a suggestion by
Kutuzov et al. (2018) to account for this fact, we aggregate the annotations
into a graded _semantic shift index_ , ranging from 0 (no shift) to 1 (shift)
depending on how many subjects spotted semantic change. The shift index is
exclusively based on the judgments of the redditors, and does not consider the
preliminary candidate selection done by us. Overall, 26 members of
r/LiverpoolFC participated in the survey, and each word received on average
8.8 judgements. Further details about the dataset are in Appendix A.

##  4 Types of Meaning Shift

We identify three main types of shift in our data via qualitative analysis of
examples with a shift index $>0.5$: metonymy, metaphor, and meme.

#### Metonymy.

In metonymic shifts, a highly salient characteristic of an entity is used to
refer to it. Among these cases are, for example, ‘highlighter’, which in
$t_{2}$ occurs in sentences like ‘we are playing with the highlighter today’,
or ‘what’s up with the hate for this kit? This is great, ten times better than
the highlighter’, used to talk about a kit in a colour similar to that of a
highlighter pen; or ‘lean’, in ‘I hope a lean comes soon!’, ‘Somebody with
speed…make a signing…Cuz I need a lean’, used to talk about hiring players due
to new hires typically leaning on a Liverpool symbol when posing for a photo
right after signing for the club. Particularly illustrative is the ‘F5’
example shown in Table 2. While ‘F5’ is initially used with its common usage
of shortcut for refreshing a page (1), it then starts to denote the act of
refreshing in order to get the latest news about the possible transfer of a
new player to LiverpoolFC (2). This use catches on and many redditors use it
to express their tension while waiting for good news (3-5),666Here the
semantic change is accompanied by a change in the part of speech, and ‘F5’
becomes a denominal verb. though not all members are aware of the new meaning
of the word (6). When the transfer is almost done, someone leaves the ‘F5
squad’ (7), and after a while, another member recalls the period in which the
word was used (8).

(1) |  Damn, after losing the F5 key on my keyboard […] | 16 Jun
---|---|---
(2) |  [he is] so close, F5 tapping is so intense right now | 18 Jun
(3) |  Don’t think about it too much, man. Just F5 | 1 Jul
(4) |  Literally 4am I slept and just woke up and thought it was f5 time | 3 Jul
(5) |  this was a happy f5 | 3 Jul
(6) |  what is an F5? | 3 Jul
(7) |  I’m leaving the f5 squad for now | 5 Jul
(8) |  I made this during the f5 madness | 6 Sep
Table 2: Examples of use of ‘F5’ with time stamps, which illustrate the speed
of the meaning shift process. All the examples are from LiverpoolFC17.

#### Metaphor.

Metaphorical shifts lead to a broadening of the original meaning of a word
through analogy. For example, in $t_{2}$ ‘shovel’ occurs in sentences such as
‘welcome aboard, here is your shovel’ or ‘you boys know how to shovel coal’:
the team is seen as a train that is running through the season, and every
supporter is asked to figuratively contribute by shoving coal into the train
boiler.

#### Meme.

Finally, memes are another prominent source of meaning shift. In this case,
fans use a word to make jokes and be sarcastic, and the new usage quickly
spreads within the community. As an example, while Liverpool was about to sign
a new player named Van Dijk, redditors started to play with the homography of
the first part of the surname with the common noun ‘van’, its plural form
‘vans’, and the shoes brand ‘Vans’: ‘Rumour has it Van Djik was wearing these
vans in the van’ or ‘How many vans could Van Dijk vear if Van Dijk could vear
vans’. Jokes of this kind are positively received in the community (‘Hahah I
love it. Anything with vans is instant karma!’) and quickly become frequent in
it.

##  5 Modeling Results and Analysis

Figure 1: Semantic shift index vs. cosine distance in the evaluation dataset
(Pearson’s $r$ = 0.49, $p<0.001$). Red horizontal ellipsis: false positives;
blue vertical ellipsis: false negatives.

The positive correlation between cosine distance and semantic shift index
(Pearson’s $r$= 0.49, $p<0.001$, see Figure 1) confirms the hypothesis that
meaning shift is mirrored by a change in context of use. However, we also find
systematic deviations.

###  5.1 False negatives

A small, but consistent group is that of words that undergo semantic shift but
are not captured by the model (blue vertical ellipsis Figure 1; shift
index>0.5, cosine distance<0.25). These are all metaphorical shifts; in
particular, cases of extended metaphor Werth (1994), where the metaphor is
developed throughout the whole text. For instance, besides the ‘shovel’
example mentioned in Section 4, we find ‘pharaoh’, the nickname of an Egyptian
player who joined Liverpool in 2017, used in sentences like ‘approved by our
new Pharaoh Tutankhamun’, or ‘our dear Egyptian Pharaoh, let’s hope he becomes
a God’. Despite the metaphoric usage, the local context of these words is
similar to the literal one, and so the model does not spot the meaning shift.
We expect this to happen in long-term shift models, too, but we are not aware
of results confirming this.

###  5.2 False positives

A larger group of problematic cases is that of words that do not undergo
semantic shift despite showing relatively large differences in context between
$t_{1}$ and $t_{2}$ (red horizontal ellipsis in Figure 1; shift index=0,
cosine distance>0.25). Manual inspection reveals that most of these “errors”
are due to a referential effect: words are used almost exclusively to refer to
a specific person or event in $t_{2}$, and so the context of use is different
from the contexts in $t_{1}$. For instance, ‘stubborn’ is almost always used
to talk about a coach who was not there in 2013 but only in 2017; ‘entourage’,
for the entourage of a particular star of the team; ‘independence’ for the
political events in Catalonia (Spain). In all these cases, the meaning of the
word stays the same, despite the change in context. In line with the
Distributional Hypothesis, the model spots the context change, but it is not
sensitive to its nature. We expect long-term shift to not be as susceptible to
referential effects like these because embeddings are aggregated over a larger
and more varied number of occurrences.

We expect that in referential cases the contexts of use will be narrower than
for words with actual semantic shift, as they are specific to one person or
event. Hence, a measure of contextual variability should help spot false
positives. To test this hypothesis, we define contextual variability as
follows: For a target word, we create a vector for each of its contexts (5
words on both sides of the target) in $t_{2}$ by averaging the embeddings of
the words occurring in it, and define variability as the average pairwise
cosine distance between context vectors.777There are alternative ways of
measuring contextual variability, but we expect them to yield the same
picture. For instance, we experimented with a different window size and
obtained the same pattern. We find that contextual variability is indeed
significantly correlated with semantic shift in our dataset (Pearson’s
$r\\!=\\!0.55$, $p\\!<\\!0.001$), while it is independent from cosine distance
(Pearson’s $r$= 0.18, $p>0.05$). These two aspects are thus complementary.
While both shift words and referential cases change context of use in $t_{2}$,
context variability captures the fact that only in referential cases words
occur in a restricted set of contexts. Figure 2 shows this effect visually.
This result can inform future models of short-term meaning shift.

Figure 2: Semantic shift index vs. context variability. Red horizontal
ellipsis: referential cases which are assigned high cosine distance values by
the model (false positives).

##  6 Conclusion

The goal of this initial study was to bring to the attention of the NLP
community short-term meaning shift, an under-studied problem in the field. Our
hope is that it will spark further research into a phenomenon which, besides
being of theoretical interest, has potential practical implications for NLP
downstream tasks concerned with user-generated language, as modeling how word
meanings rapidly change in communities would allow a better understanding of
what their members say. Future research should experiment with other datasets
(reddits from other domains, other online communities) and also alternative
models that address the challenges described here.

## Acknowledgements

The research carried out by the Amsterdam section of the team was partially
funded by the Netherlands Organisation for Scientific Research (NWO) under
VIDI grant no. 276-89-008, Asymmetry in Conversation. This project has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement No
715154), and from the Spanish Ramón y Cajal programme (grant RYC-2015-18907).
This paper reflects the authors’ view only, and the EU is not responsible for
any use that may be made of the information it contains.

## References

  * Artstein and Poesio (2008) Ron Artstein and Massimo Poesio. 2008.  Inter-coder agreement for computational linguistics.  _Computational Linguistics_ , 34(4):555–596.
  * Azarbonyad et al. (2017) Hosein Azarbonyad, Mostafa Dehghani, Kaspar Beelen, Alexandra Arkut, Maarten Marx, and Jaap Kamps. 2017.  Words are malleable: Computing semantic shifts in political and media discourse.  In _Proceedings of the 2017 ACM on Conference on Information and Knowledge Management_ , pages 1509–1518. ACM.
  * Brennan and Clark (1996) Susan E. Brennan and Herbert H. Clark. 1996.  Conceptual pacts and lexical choice in conversation.  _Journal of Experimental Psychology: Learning, Memory, and Cognition_ , 22(6):1482.
  * Clark (1996) Herbert H. Clark. 1996.  _Using language_.  Cambridge University Press.
  * Danescu-Niculescu-Mizil et al. (2013) Cristian Danescu-Niculescu-Mizil, Robert West, Dan Jurafsky, Jure Leskovec, and Christopher Potts. 2013.  No country for old members: User lifecycle and linguistic change in online communities.  In _Proceedings of the 22nd international conference on World Wide Web_ , pages 307–318. ACM.
  * Del Tredici and Fernández (2017) Marco Del Tredici and Raquel Fernández. 2017.  Semantic variation in online communities of practice.  In _IWCS 2017-12th International Conference on Computational Semantics-Long papers_.
  * Del Tredici and Fernández (2018) Marco Del Tredici and Raquel Fernández. 2018.  The road to success: Assessing the fate of linguistic innovations in online communities.  In _COLING 2018-27th International Conference on Computational Linguistics-Long papers_.
  * Del Tredici et al. (2016) Marco Del Tredici, Malvina Nissim, and Andrea Zaninello. 2016.  Tracing metaphors in time through self-distance in vector spaces.  In _CLiC it 2016 - Third Italian Conference on Computational Linguistics_.
  * Eckert and McConnell-Ginet (1992) Penelope Eckert and Sally McConnell-Ginet. 1992.  Communities of practice: Where language, gender, and power all live.  In _Locating Power, Proceedings of the 1992 Berkeley Women and Language Conference_ , pages 89–99.
  * Frermann and Lapata (2016) Lea Frermann and Mirella Lapata. 2016.  A bayesian model of diachronic meaning change.  _Transactions of the Association for Computational Linguistics_ , 4:31–45.
  * Gulordava and Baroni (2011) Kristina Gulordava and Marco Baroni. 2011.  A distributional similarity approach to the detection of semantic change in the google books ngram corpus.  In _Proceedings of the GEMS 2011 Workshop on GEometrical Models of Natural Language Semantics_ , pages 67–71. Association for Computational Linguistics.
  * Hamilton et al. (2016) William L Hamilton, Jure Leskovec, and Dan Jurafsky. 2016.  Diachronic word embeddings reveal statistical laws of semantic change.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , volume 1, pages 1489–1501.
  * Harris (1954) Zellig S Harris. 1954.  Distributional structure.  _Word_ , 10(2-3):146–162.
  * Hasan (2009) Ruqaiya Hasan. 2009.  _Semantic variation: Meaning in society and in sociolinguistics_.  Equinox London.
  * Jatowt and Duh (2014) Adam Jatowt and Kevin Duh. 2014.  A framework for analyzing semantic change of words across time.  In _Proceedings of the 14th ACM/IEEE-CS Joint Conference on Digital Libraries_ , pages 229–238. IEEE Press.
  * Kim et al. (2014) Yoon Kim, Yi-I Chiu, Kentaro Hanaki, Darshan Hegde, and Slav Petrov. 2014.  Temporal analysis of language through neural language models.  In _Proceedings of the ACL 2014 Workshop on Language Technologies and Computational Social Science_ , pages 61–65. Association for Computational Linguistics.
  * Kulkarni et al. (2015) Vivek Kulkarni, Rami Al-Rfou, Bryan Perozzi, and Steven Skiena. 2015.  Statistically significant detection of linguistic change.  In _Proceedings of the 24th International Conference on World Wide Web_ , pages 625–635. ACM.
  * Kutuzov et al. (2018) Andrey Kutuzov, Lilja Øvrelid, Terrence Szymanski, and Erik Velldal. 2018.  Diachronic word embeddings and semantic shifts: a survey.  In _Proceedings of the 27th International Conference on Computational Linguistics (COLING)_ , pages 1384–1397. Association for Computational Linguistics.
  * Levy et al. (2015) Omer Levy, Yoav Goldberg, and Ido Dagan. 2015.  Improving distributional similarity with lessons learned from word embeddings.  _Transactions of the Association for Computational Linguistics_ , 3:211–225.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013.  Distributed representations of words and phrases and their compositionality.  In _Advances in neural information processing systems_ , pages 3111–3119.
  * Phillips et al. (2017) Lawrence Phillips, Kyle Shaffer, Dustin Arendt, Nathan Hodas, and Svitlana Volkova. 2017.  Intrinsic and extrinsic evaluation of spatiotemporal text representations in twitter streams.  In _Proceedings of the 2nd Workshop on Representation Learning for NLP_ , pages 201–210.
  * Rotabi et al. (2017) Rahmtin Rotabi, Cristian Danescu-Niculescu-Mizil, and Jon Kleinberg. 2017.  Competition and selection among conventions.  In _Proceedings of the 26th International Conference on World Wide Web_ , pages 1361–1370. International World Wide Web Conferences Steering Committee.
  * Sagi et al. (2011) Eyal Sagi, Stefan Kaufmann, and Brady Clark. 2011.  Tracing semantic change with latent semantic analysis.  _Current methods in historical semantics_ , pages 161–183.
  * Stewart and Eisenstein (2018) Ian Stewart and Jacob Eisenstein. 2018.  Making” fetch” happen: The influence of social and linguistic context on the success of lexical innovations.  _Proceedings of Conference on Empirical Methods in Natural Language Processing (EMNLP)_.
  * Szymanski (2017) Terrence Szymanski. 2017.  Temporal word analogies: Identifying lexical replacement with diachronic word embeddings.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , volume 2, pages 448–453.
  * Tang (2018) Xuri Tang. 2018.  A state-of-the-art of semantic change computation.  _Natural Language Engineering_ , 24(5):649––676.
  * Wenger (1998) Etienne Wenger. 1998.  _Communities of practice: Learning, meaning, and identity_.  Cambridge University Press.
  * Werth (1994) Paul Werth. 1994.  Extended metaphor—a text-world account.  _Language and literature_ , 3(2):79–103.
  * Wijaya and Yeniterzi (2011) Derry Tanti Wijaya and Reyyan Yeniterzi. 2011.  Understanding semantic change of words over centuries.  In _Proceedings of the 2011 international workshop on DETecting and Exploiting Cultural diversiTy on the social web_ , pages 35–40. ACM.

##  Appendix A Further Details on Evaluation Dataset

For our experiment, we considered content words only, which we identified by
using the external list of common words available at
https://www.wordfrequency.info/free.asp.

Three words were discarded from the initial list after analysis of the
redditor data: ‘discord’ and ‘owls’ due to the homonymy with proper names not
detected during survey’s implementation; ‘tracking’ because the chosen
examples clearly mislead the judgements of the redditors.

As detailed in Section 3.3, 26 members of r/LiverpoolFC participated in the
survey, and each word received on average 8.8 judgements. We computed inter-
annotator agreement as Krippendorff’s alpha, and obtained $\alpha$ = 0.58, a
relatively low value but common in semantic tasks Artstein and Poesio (2008).

The results of the annotation validate our initial word sampling procedure:

  * •

the words that present a significant increase in frequency and were annotated
as meaning shift by us received an average shift annotation of 0.72 ($\pm$
0.15);

  * •

the words that present a significant increase in frequency but that were _not_
annotated as meaning shift by us received an average shift annotation of 0.15
($\pm$ 0.16);

  * •

the words that keep a constant frequency between $t_{1}$ and $t_{2}$, and we
don’t consider examples of meaning shift, got 0.07 ($\pm$ 0.12).