\documentclass[cmbright,fleqn,referee]{envauth}
\usepackage{multirow}
\usepackage[FIGTOPCAP]{subfigure}
%%%%% AUTHORS - PLACE YOUR OWN MACROS HERE %%%%%
%\newcommand{\tr}{\top} % Transpose
\newcommand{\chp}[1]{\widehat{#1}} % Hat
\newcommand{\argmax}[1]{\underset{#1}{\operatorname{argmax}}}
\newcommand{\bE}{\textrm{E}}
\newcommand{\bmu}{\mbox{\boldmath $\mu$}}
\newcommand{\bphi}{\mbox{\boldmath $\phi$}}
\newcommand{\bPhi}{\mbox{\boldmath $\Phi$}}
\newcommand{\bnu}{\mbox{\boldmath $\nu$}}
\newcommand{\bSigma}{\mbox{\boldmath $\Sigma$}}
\newcommand{\bepsilon}{\mbox{\boldmath $\epsilon$}}
\newcommand{\bLambda}{\mbox{\boldmath $\Lambda$}}
\newcommand{\bbeta}{\mbox{\boldmath $\beta$}}
\newcommand{\btheta}{\mbox{\boldmath $\theta$}}
\newcommand{\bTheta}{\mbox{\boldmath $\Theta$}}
\newcommand{\bxi}{\mbox{\boldmath $\xi$}}
\newcommand{\bvarphi}{\mbox{\boldmath $\varphi$}}
\newcommand{\bnabla}{\mbox{\boldmath $\nabla$}}
\newcommand{\bDelta}{\mbox{\boldmath $\Delta$}}
\newcommand{\bdelta}{\mbox{\boldmath $\delta$}}
\newcommand{\bvar}{\mbox{\boldmath $\varepsilon$}}
\newcommand{\bzeta}{\mbox{\boldmath $\zeta$}}
\newcommand{\NI}{\textrm{NI}}
\newcommand{\N}{\textrm{N}}
\newcommand{\SNI}{\textrm{SNI}}
\newcommand{\SN}{\textrm{SN}}
\newcommand{\ST}{\textrm{ST}}
\newcommand{\SSL}{\textrm{SSL}}
\newcommand{\SCN}{\textrm{SCN}}
\newcommand{\SMN}{\textrm{SMN}}
\newcommand{\A}{\mathbf{A}}
\newcommand{\B}{\mathbf{B}}
\newcommand{\bV}{\mathbf{V}}
\newcommand{\bPsi}{\mbox{\boldmath $\Psi$}}
\newcommand{\bpsi}{\mbox{\boldmath $\psi$}}
\newcommand{\btau}{\mbox{\boldmath $\tau$}}
\newcommand{\bpi}{\mbox{\boldmath $\pi$}}
\newcommand{\bOmega}{\mbox{\boldmath $\Omega$}}
\newcommand{\bomega}{\mbox{\boldmath $\omega$}}
\newcommand{\bupsilon}{\mbox{\boldmath $\upsilon$}}
\newcommand{\bgamma}{\mbox{\boldmath $\gamma$}}
\newcommand{\bGamma}{\mbox{\boldmath $\Gamma$}}
\newcommand{\neta}{\mbox{\boldmath $\eta$}}
\newcommand{\be}{\mathbf{b}}
\newcommand{\bF}{\mathbf{F}}
\newcommand{\bQ}{\mathbf{Q}}
\newcommand{\bP}{\mathbf{P}}
\newcommand{\bC}{\mathbf{C}}
\newcommand{\yp}{\mathbf{y}}
%\newcommand{\xp}{\mathbf{x}}
\newcommand{\x}{\mathbf{x}}
\newcommand{\y}{\mathbf{y}}
\newcommand{\Y}{\mathbf{Y}}
\newcommand{\vari}{\phi_b\textbf{u}_j\textbf{u}'_j+\mathbb{I}_{n_j}\phi_e}
\newcommand{\bD}{\mathbf{D}}
\newcommand{\C}{\mathbf{C}}
\newcommand{\blambda}{\mbox{\boldmath $\lambda$}}
\newcommand{\brho}{\mbox{\boldmath $\rho$}}
\newcommand{\z}{\mathbf{Z}}
\newcommand{\zp}{\mathbf{z}}
\newcommand{\Z}{\mathbf{Z}}
\newcommand{\R}{\mathbf{R}}
\newcommand{\X}{\mathbf{X}}
\newcommand{\I}{\mathbf{I}}
\newcommand{\bW}{\mathbf{W}}
\newcommand{\bH}{\mathbf{H}}
\newcommand{\E}{\textrm{E}}
\newcommand{\ba}{\mathbf{b}}
\newcommand{\bd}{\mathbf{d}}
\newcommand{\by}{{\bf y}}
\newcommand{\sa}{\mathbf{s}}
\newcommand{\bu}{\mathbf{u}}
\newcommand{\bS}{\mathbf{S}}
\newcommand {\mbf}[1] {\mbox {\boldmath$#1$} }
\newcommand{\ii}{i=1,\ldots,n}
\newcommand{\jj}{j=1,\ldots,p}
%\newcommand{\by}{{\bf Y}_j}
\newcommand{\balpha}{\mbox{\boldmath $\alpha$}}
\newcommand{\qc}{\stackrel {{\rm q.c}}{\longrightarrow}}
\newcommand{\iid}{\stackrel {{\rm iid.}}{\sim}}
\newcommand{\ind}{\stackrel {{\rm ind.}}{\sim}}
\newcommand{\app}{\stackrel {{\rm .}}{\sim}}
\newcommand{\dd}{\stackrel {{\rm d}}{\longrightarrow}}
\newcommand{\sumas}{\sum^n_{i=1}}
\newcommand{\sumasp}{\sum^{p}_{j=1}}
%\newcommand{\tr}{\textrm{tr}}
%\newtheorem{theorem}{Theorem}[section]
\renewcommand{\baselinestretch}{1.3}
\newcommand{\s}{\mathbf{s}}
\newcommand{\PR}[1]{\ensuremath{\left[#1\right]}}
\newcommand{\PC}[1]{\ensuremath{\left(#1\right)}}
\newcommand{\chav}[1]{\ensuremath{\left\{#1\right\}}}
\newcommand{\V}{\mathbf{V}}
\def\bSig\mathbf{\Sigma}
\newcommand{\VS}{V\&S}
\newcommand{\tr}{\mbox{tr}}
\received{XX Month 201X}
\revised{XX Month 201X}
\accepted{XX Month 201X}
\newtheorem{theorem}{Theorem}
\usepackage{natbib,upgreek}
%%%%%%%
\usepackage{graphicx}
\usepackage{epstopdf}
%%%%%%%
%\usepackage[nolists]{endfloat}
\runninghead{Lachos et al.}{Diagnostics in Spatial Models}
\begin{document}
\title{Influence Diagnostics in Spatial Models with Censored Response}
\author{Victor H. Lachos\affil{a}\corrauth, Larissa A. Matos\affil{a}, Thais S. Barbosa \affil{a}, \\ Aldo M. Garay\affil{b} and Dipak K. Dey\affil{c}}
\corraddr{V\'{\i}ctor Hugo Lachos D\'{a}vila, Departamento de Estat\'{\i}stica, IMECC, Universidade Estadual de Campinas, CEP 13083-859, Campinas, S\~{a}o Paulo, Brazil. E-mail: hlachos@ime.unicamp.br}
\address{\affilnum{a}Department of Statistics, Campinas State University, Campinas, S\~{a}o Paulo, Brazil\\
\affilnum{b}Department of Statistics, Federal University of Pernambuco, Recife, Brazil\\
\affilnum{c}Department of Statistics, University of Connecticut, Storrs, CT, USA
}
\begin{abstract}
Environmental data are often spatially correlated and sometimes include observations below or above detection limits (i.e., censored values reported as less or more than a level of detection). Existing research mainly concentrate on parameter estimation using Gibbs sampling, and most researches conducted from a frequentist perspective in spatial censored models are elusive. In this paper, we propose an exact estimation procedure to obtain the maximum likelihood estimates of the fixed effects and variance components, using a stochastic approximation of the EM (SAEM) algorithm \citep{delyon1999convergence}. This approach permits estimation of the parameters of spatial linear models when censoring is present in an easy and fast way. As a byproduct, predictions of unobservable values of the response variable are possible. Motivated by this algorithm, we develop local and global influence measures on the basis of the conditional expectation of the complete-data log-likelihood function which eliminates the complexity associated with the approach of \cite{cook1977detection,cook86} for spatial censored models. Some useful perturbation schemes are discussed. The newly developed method is illustrated using data from a dioxin contaminated site in Missouri that contain left-censored data as well as a dataset related to depths of a geological horizon that contains both left- and right-censored observations. In addition, a simulation study is presented that, explores the accuracy of the proposed measures in detecting influential observations under different perturbation schemes. The methodology addressed in this paper is implemented in the R package {\it CensSpatial}.
\end{abstract}
\keywords{Censored data; Geostatistical data; SAEM Algorithm; Influential observations.}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}\label{s:intro}
Spatial data are common in ecology, environmental health, mining, hydrology and epidemiology, where sampling units are geographical areas or spatially located individuals. Analysis of spatial data is challenged by the spatial correlation among the observations such as the conditional autoregressive (CAR) structure or the Mat\'{e}rn correlation structure. An additional complication is that spatial data are subject to upper or lower detection limits below and above which they are not quantifiable. For example, environmental (spatial) monitoring of different variables often involves left-censored observations falling below the minimum limit of detection (LOD) of the instruments used to quantify them. The proportion of censored data in these studies may not be small, so the use of crude/ad hoc methods, such as substituting a threshold value or some arbitrary point like a midpoint between zero and cutoff for detection, might lead to biased estimates of fixed effects and variance components \citep{de2005bayesian,fridley2006data}.
As an alternative to crude imputation methods \cite{militino1999analyzing} develop an EM algorithm for maximum likelihood (ML) estimation in censored spatial data. However, this approach suffers from several drawbacks that restrict its applicability. For instance, \cite{de2005bayesian} notes that this ML approach does not provide a means to estimate the correlation structure in the data and hence assume it is known. In addition, when estimating variables at non-sampled locations the observed and imputed data are not differentiated, thus underestimating the predictive uncertainty. Due to the complexity of the likelihood-based methods, which involve computationally intractable integrals, \cite{de2005bayesian} and \cite{fridley2006data} adopt a Bayesian approach to inference and prediction for spatially correlated censored observations. In both papers, data augmentation and Markov chain Monte Carlo (MCMC) algorithms are utilized. In this paper, we first propose a stochastic version of the EM algorithm for ML estimation, the so-called SAEM algorithm proposed by \cite{delyon1999convergence}. Then, the diagnostic measures for assessing the local influence in spatial censored linear (SCL) models are developed and presented. In the framework of spatial models, \cite{jank2006implementing} shows that the computational effort of SAEM is much smaller and reaches convergence in just a fraction of the simulation size when compared to the Monte Carlo EM (MCEM) algorithm. This is due to the memory effect contained in the SAEM method, in which the previous simulations are considered in the computation of the posterior ones.
The study of influential values is an important and a key step in data analysis subsequent to parameter estimation. This can be carried out by conducting {an influence analysis to detect influential observations. There are two primary approaches for detecting influential observations. The first approach is the case-deletion approach \citep{cook1977detection} and is an intuitively appealing method (see also Cook and Weisberg, 1982). Deletion diagnostics such as Cook's distance or the likelihood distance have been applied to many statistical models. The second approach, which is a general statistical technique used to assess the stability of the estimation outputs with respect to the model inputs, is the local influence approach of \cite{cook86}}. Following the pioneering work of \cite{cook86}, this method has received considerable attention recently in the statistical literature of spatial models; see, for example, \cite{assumpccao2014analysis} and \cite{de2014influence}.
Although several diagnostic studies on spatial models have appeared in the literature, to the best of our knowledge, no study has been made conducted of influence diagnostics for censored spatial data or of the local influence analysis. The main difficulty is due to the fact that the observed log-likelihood function of censored spatial models involves intractable integrals (for instance, the pdfs of truncated multinormal distributions), making the direct application of Cook's approach \citep{cook86} to this model very difficult if not impossible, since the measures involve the first and second partial derivatives of this function. \cite{ZhuLee2001} develop an approach for performing local influence analysis for general statistical models with missing data, based on the {\it Q}-displacement function which is closely related to the conditional expectation of the complete-data log-likelihood in the E-step of the EM algorithm. {{This approach produces results very similar to those obtained from Cook's method. Moreover, the case-deletion can be studied by {\it Q}-displacement function following the approach of \cite{zhu2001case}. So, we develop here methods to obtain case-deletion measures and local influence measures by using the method of \cite{zhu2001case} \citep[see also][]{ZhuLee2001, LeeXu04} in the context of SCL models.
The rest of the paper is organized as follows. Section \ref{model:est} gives a brief description of the spatial linear model, including an outline of the SAEM algorithm. Section \ref{CensSpat} proposes the SCL model and shows how to get the ML estimates through the SAEM algorithm. In Section \ref{a:diagnostic}, we provide a brief sketch of the local influence approach for models with incomplete data, and also develop a method pertinent to the SCL model. Three different perturbation schemes are considered. The method is illustrated in Section \ref{realdata} with the analysis of a dataset from depths of a geological horizon containing left- and right- censored data, as well as the analysis of a dataset from a dioxin contaminated site in Missouri containing left censored data. Section \ref{simulation} presents simulation studies. Finally, Section \ref{concl} concludes with a short discussion of issues raised by our study and some possible directions for a future research.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Preliminaries}\label{model:est}
\subsection{\textbf{The Spatial linear model}}\label{subSAEM}
As presented by \cite{de2014influence}, we consider a Gaussian stochastic process $\{Z(\sa), \sa\in \mathbb{D}\}$, where $\mathbb{D}$ is a subset of $\mathbb{R}^d$ , the $d$-dimensional Euclidean space. It supposes that data $Z(\sa_1), . . . , Z(\sa_n)$ of this process are observed at known sites (locations) $\sa_i$, for $\ii$, where $\sa_i$ is a $d$-dimensional vector of spatial site coordinates, and are generated from the model:
\begin{equation}
Z(\sa_i)=\mu(\sa_i)+\epsilon(\sa_i),
\end{equation}
where both the deterministic term $\mu(\sa_i)$ and the stochastic term $\epsilon(\sa_i)$ can depend on the spatial location at which $Z(\sa_i)$ is observed. We assume that the stochastic errors have zero mean, $E\{\epsilon(\sa_i)\}=0$, and that variation between spatial points is determined by a covariance function $C(\sa_i,\sa_j)=\textrm{Cov}\{\epsilon(\sa_i), \epsilon(\sa_j)\}$. Suppose that for some known functions of $\sa_i$, $x_1(\sa_i),\ldots, x_p(\sa_i)$, the mean of the stochastic process is
\begin{equation}\label{mu}
\mu(\sa_i)=\sumasp x_j(\sa_i)\beta_j,
\end{equation}
where $\bbeta=(\beta_1,\ldots,\beta_p)^{\top}$ are unknown parameters to be estimated. In addition, each family of covariance functions $C(\sa_i,\sa_j)$, is fully specified by a $q$-dimensional parameter vector $\bphi=(\phi_1,\ldots,\phi_q)^{\top}$.
Considering $\mu(\sa_i)=\x_i^{\top}\bbeta$ and $Z_i=\x_i^{\top}\bbeta+\epsilon_i$, for $\ii$, the spatial linear model is given by
\begin{equation}\label{matrixSpatial}
\Z=\X\bbeta+\bepsilon,
\end{equation}
where $\Z=(Z_1,\ldots,Z_n)^{\top}~$ with $Z_i=Z(\sa_i),~$ $\X$ is the $n\times p$ matrix (full rank) with $i-$th row $\x^{\top}_i=(x_{i1},\ldots,x_{ip})^{\top}$ and $x_{ij}=x_j(\sa_i)$ for $\ii$ and $\jj$.
Note that $\bepsilon=(\epsilon_1,\ldots,\epsilon_n)^{\top}$ with $\epsilon_i=\epsilon(\sa_i)$, where $E\{\bepsilon\}=\textbf{0}$ and the scale matrix of $\bepsilon$ is the non singular matrix $\bSigma=[C(\sa_i,\sa_j)]=\tau^2\mathbf{I}_n+\sigma^2\R(\rho)$. Also, $\tau^2$ can be viewed as a measurement error variance or a nugget effect, $\sigma^2$ is defined as the sill, $\R=\R(\rho)=[r_{ij}],$ is an $n\times n$ symmetric matrix with diagonal elements $r_{ii}=1$, for $\ii$ and $\rho$ is a function of the range of the model. In general, $\R$ depends on the Euclidean distance $d_{ij} = ||\sa_i - \sa_j||$ between the points $\sa_i$ and $\sa_j$. This
parametric form occurs in several isotropic processes, as presented by Table \ref{Corr_Isot}. It is important to note that the spherical family lacks flexibility in comparison with the two-parameter Mat\'{e}rn class. Also, $\R(\rho) $ is only once differentiable at $d_{ij}=\rho$, which causes technical difficulties with maximum likelihood estimation.
\begin{table}[h!]
\centering \large \caption{Different types of correlation function
$\R$.} \label{Corr_Isot}
\begin{center}
\scriptsize{
%\footnotesize{
%\large{
\begin{tabular*}{\textwidth}{l @{\extracolsep{\fill}} |l|l}
\toprule
%\scriptsize
\textbf{Parametric forms } &$\R(\rho)=\R(\rho,d_{ij})$ & \scriptsize \textbf{Description} \\ \midrule
\multirow{3}{*}{\textbf{(1)} Mat\'{e}rn Family } & \scriptsize \multirow{3}{*}{$\left\{\begin{array}{ll}
\displaystyle \frac{1}{2^{\kappa-1}\Gamma{(\kappa)}}\left(\frac{d_{ij}}{\rho}\right)^{\kappa}K_{\kappa}(d_{ij}/\rho),& d_{ij}>0, \\
1, & d_{ij}=0, \end{array}\right.$} & $\rho>0$ and $\kappa > 0$ fixed. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \\
& & $K_{\kappa}(u)=\frac{1}{2}\int_{0}^{\infty}x^{\kappa-1}e^{-\frac{1}{2}u(x+x^{-1})}dx$ \\
& & \\ \hline
\multirow{3}{*}{\textbf{(1.1)} Gaussian covariance} & \scriptsize \multirow{3}{*}{$\left\{\begin{array}{ll}
\exp\left\{-\left(\displaystyle\frac{d_{ij}}{\rho}\right)^2\right\},&d_{ij}>0, \\
1,\,\,\,& d_{ij}=0.
\end{array}\right.$} & \\
& & A special case of Mat\'{e}rn family when $\kappa\rightarrow\infty$ \\
& & \\ \hline
\multirow{3}{*}{\textbf{(1.2)} Exponential covariance } &\scriptsize \multirow{3}{*}{$\left\{\begin{array}{ll}
\exp\left\{-\left(\displaystyle\frac{d_{ij}}{\rho}\right)\right\},&d_{ij}>0, \\
1,&d_{ij}=0. \end{array}\right.$} & \\
& & A special case of Mat\'{e}rn family when $\kappa=1/2$\\
& & \\ \hline
\multirow{3}{*}{\textbf{(2)} Spherical Family} & \scriptsize \multirow{3}{*}{$\left\{\begin{array}{ll}
1-\displaystyle\frac{3}{2}(d_{ij}/\rho)+\displaystyle\frac{1}{2}(d_{ij}/\rho)^3,& 0\leq d_{ij} \leq \rho, \\
0,&d_{ij}>\rho.
\end{array}\right.$} & $\rho > 0$ is a single parameter. \\
& & \\
& & $\R(\rho) = 0$ for sufficiently large $d_{ij}$, namely $ d_{ij} > \rho$. \\
\bottomrule
\end{tabular*}
}
\end{center}
\end{table}
%{\color{blue}
\subsection{The EM and SAEM algorithms}\label{secSAEM}
In models with missing and censored data, the EM algorithm
\citep{Dempster77} has become established as the most popular tool
to obtain the ML estimates of model parameters. Define
$\mathbf{Z}_{\mathrm{com}}=(\mathbf{Z}_{\mathrm{m}},\mathbf{Z}_{\mathrm{obs}})$,
where $\mathbf{Z}_{\mathrm{m}}$ denotes the missing data and
$\mathbf{Z}_{\mathrm{obs}}$ the observed data. This iterative
algorithm maximizes the complete log-likelihood function
$\ell_{c}(\mathbf{\btheta};\,\mathbf{Z}_{\mathrm{com}})$ at each
step, converging quickly to a stationary point of the observed
likelihood $(\ell(\mathbf{\btheta};\,\mathbf{Z}_{\mathrm{obs}}))$
under mild regularity conditions
\citep{wu1983convergence,vaida2005parameter}.
The EM algorithm
proceeds in two simple steps, \textbf{\emph{E-Step:}} Replace the observed likelihood by
the complete likelihood and compute its conditional expectation
$\mathit{Q}(\mathbf{\btheta}|\widehat{\mathbf{\btheta}}^{(k)})=\mathit{\mathrm{E}}\left\{
\ell_{c}(\mathbf{\btheta};\,\mathbf{Z}_{\mathrm{com}})|\widehat{\mathbf{\btheta}}^{(k)},\mathbf{Z_{\mathrm{obs}}}\right\}
,$ where $\widehat{\mathbf{\btheta}}^{(k)}$ is the estimate of
$\btheta$ at the \emph{k}-th iteration;
\textbf{\emph{M-Step:}} Maximize
$\mathit{Q}(\mathbf{\theta}|\widehat{\mathbf{\btheta}}^{(k)})$ with
respect to $\btheta$ to obtain $\widehat{\mathbf{\btheta}}^{(k+1)}$.
However, in some applications of the EM algorithm, the E-step cannot
be obtained analytically, so it has to be calculated using simulations.
\cite{wei1990monte} proposed the Monte Carlo EM (MCEM) algorithm in
which the E-step is replaced by a Monte Carlo approximation based on
a large number of independent simulations of the missing data. But this
simple solution is computationally expensive, given the need
to generate a large number of independent simulations of the missing
data for a good approximation. Thus, in order to reduce the amount
of required simulations compared to the MCEM algorithm, the SAEM
algorithm proposed by \cite{delyon1999convergence} replaces the
E-step of the EM algorithm by a stochastic approximation procedure,
while the maximization step remains unchanged. Besides having good
theoretical properties, the SAEM estimates the population parameters
accurately, converging to the global maxima of the ML estimates
under quite general conditions
\citep{allassonniere2010construction,delyon1999convergence,kuhn2004coupling}.
At each iteration, the SAEM algorithm successively simulates missing
data with the conditional distribution, and
updates the unknown parameters of the model. Thus, at iteration \emph{k}, the SAEM proceeds as follows \\
\textbf{\textit{E-Step:}} (1)
Simulation: Draw $(\mathbf{q}^{(\ell,k)}),\, \ell=1,\ldots,m$ from the conditional distribution $f(\mathbf{q}|\mathbf{\mathbf{\btheta}}^{(k-1)},\mathbf{Z}_{i})$; and
(2) Stochastic Approximation: Update the $\mathit{Q}(\mathbf{\theta}|\widehat{\mathbf{\theta}}^{(k)})$ function as
\begin{eqnarray}\label{saem0}
\mathit{Q}(\mathbf{\btheta}|\widehat{\mathbf{\btheta}}^{(k)})\approx\mathit{Q}(\mathbf{\btheta}|\widehat{\mathbf{\btheta}}^{(k-1)})+\delta_{k}\left[\frac{1}{m}\sum_{\ell=1}^{m}\ell_{c}(\mathbf{\theta}|\mathbf{Z}_{\mathrm{obs}},\mathbf{q}^{(\ell,k)})\mathit{-}\mathit{Q}(\mathbf{\btheta}|\widehat{\mathbf{\btheta}}^{(k-1)})\right],
\end{eqnarray}
where $\delta_{k}$ is a smoothness parameter.
\noindent\textbf{\emph{M-Step}} (Maximization): Update $\widehat{\mathbf{\btheta}}^{(k)}$ as $\widehat{\mathbf{\btheta}}^{(k+1)}=\argmax{\btheta}\mathit{Q}(\mathbf{\btheta}|\widehat{\mathbf{\btheta}}^{(k)})$.
An important step on SAEM algorithm is to define its parameters. As mentioned by \cite{kuhn2004coupling},
$\delta_{k}$ is a smoothness parameter, i.e., a decreasing sequence of positive numbers such that $\sum_{k=1}^{\infty}\delta_{k}=\infty$ and $\sum_{k=1}^{\infty}\delta_{k}^{2}<\infty$. Note that, for the SAEM algorithm, the E-Step coincides with the MCEM algorithm, but only a small number of simulations $m$ (suggested to be $m\leq20$) is necessary. This is possible because unlike the traditional EM algorithm and its variants, the SAEM algorithm uses not only the current simulation of the missing data at iteration \emph{k} denoted by $\mathbf{q}^{\scriptscriptstyle(\ell,k)},\, \ell=1,\ldots,m$, but some or all previous simulations, where this ``memory'' property is set by the smoothing parameter $\delta_{k}$.
Note that in Equation (\ref{saem0}), if the smoothing parameter $\delta_{k}$ is equal to 1 for all $k$, the SAEM algorithm will have ``no memory'', and will be equivalent to the MCEM algorithm. The SAEM with no memory will converge quickly (convergence in distribution) to a solution neighborhood, but the algorithm with memory will converge slowly (almost sure convergence) to the ML solution. As proposed by \cite{Galarza2015}, we considered the following choice of the smoothing parameter
\[ \delta_{k}=\begin{cases}
1, & \quad\mathrm{for}\quad\quad\quad\,\,1\leq k\leq cW,\\
\frac{1}{k-cW}, & \quad\mathrm{for}\quad cW+1\leq k\leq W,
\end{cases}
\]
where $W$ is the maximum number of iterations, and $c$ a cutoff point ($0\leq c\leq1$) which determines the percentage of initial iterations with no memory. For example, if $c=0$, the algorithm will have memory for all iterations, and hence will converge slowly to the ML estimates. If $c=1$, the algorithm will have no memory, and so will converge quickly to a solution neighborhood. For the first case, $W$ would need to be large in order to achieve the ML estimates. For the second, the algorithm will put out a Markov Chain where after applying a \emph{burn in} and \emph{thinning}, the mean of the chain observations can be a reasonably estimated.
A number between 0 and 1 ($0\zeta_{r+1}=\ldots=0$ and orthonormal eigenvectors $\{\bvar_k, k=1,\ldots,g \}$, \cite{ZhuLee2001} propose to inspect all eigenvectors corresponding to nonzero eigenvalues for capturing more information. Following the work of \cite{ZhuLee2001}, we consider the following aggregated contribution vector of all eigenvectors that correspond to nonzero eigenvalues. Let $\tilde{\zeta}_k=\zeta_k/(\zeta_1+\ldots+\zeta_r),$ $\bvar^2_k=(\bvar_{k1}^2,\ldots,\bvar_{kg}^2)^{\top}$ and $M(0)=\sum_{k=1}^r\tilde{\zeta}_k\bvar_k^2.$ The $l-th$ component of $M(0)$, $M(0)_l$, is equal to $\sum_{k=1}^r\tilde{\zeta}_k\bvar_{kl}^2$. The assessment of influential cases is based on visual inspection of $\{M(0)_l,l=1,\ldots,g\}$ plotted against the index $l$. The $l-th$ case may be regarded as influential if $M(0)_l$ is larger than the benchmark value.}}
{{The inconvenience in using the normal curvature is deciding about the influence of the observations, since $C_{f_Q,\textbf{d}}(\btheta)$ may assume any value and it is not invariant under a uniform change of scale. Based on the work of \cite{poon1999conformal} in using a conformal normal curvature, \cite{ZhuLee2001} considered the conformal normal curvature {{$B_{f_Q,\textbf{d}}(\btheta)=C_{f_Q,\textbf{d}}(\btheta)/tr[-2\ddot{Q}_{\bomega_0}],$ whose computation is quite simple and also has the property that $0 \leq B_{f_Q,\textbf{d}}(\btheta) \leq 1$.}} Let $\textbf{d}_l$ be a basic perturbation vector with $l$th entry as $1$ and all other entries as 0. \cite{ZhuLee2001} then showed that for all $l$, {{$M(0)_l=B_{f_Q,\textbf{d}_l}$}}. Therefore, we can obtain $M(0)_l$ via {{ $B_{f_Q,\textbf{d}_l}$}}.}}
{{So far, there is no general rule to judge how large the influence of a specific case in the data is. Let $\overline{M}(0)$ and $SM(0)$ denote, respectively, the mean and standard error of $\{M(0)_l: l=1,\ldots,g\},$ where $\overline{M(0)}=1/g$. \cite{poon1999conformal} propose using $2\overline{M}(0)$ as a benchmark for $M(0)$. For instance, \cite{ZhuLee2001} propose using $\overline{M}(0)+ 2 SM(0)$ as a benchmark to take into account the variance of $\{M(0)_l: l=1,\ldots,g\}$ as well. According to \cite{LeeXu04}, the exact choice of the function of $\overline{M}(0)$ as the benchmark is subjective. \cite{LeeXu04} also propose using $\overline{M}(0)+ c^{*} SM(0)$, where $c^{*}$ is a selected constant. The choice of $c^{*}$ is subjective.
In this paper, we will consider $c^{*} = 3$; following \cite{matos2013influence}.
%{\color{blue}
%}
\subsection{\textbf{Perturbation scheme}}\label{s:PerScheme}
In this section, we will evaluate the $\bDelta$ matrix under the following perturbation schemes for {SCL} models: {\it Perturbation of response variables} carried out on the response values, which may indicate observations with large influence on their own predicted values (in our case, the response variables are $\bV's$); {\it Scale perturbation} performed on the scale matrix $\bSigma=\tau^2\mathbf{I}_{n}+\sigma^2R(\rho)$, which may reveal individuals that are most influential on the scale structure and consequently on the $\balpha$ estimate and finally \emph{perturbation of explanatory variables}.
{Let $\btheta= (\bbeta^{\top}, \balpha^{\top})^{\top}$, where $\balpha=(\sigma^2,\tau^2,\rho)^{\top}$. Given that, the matrix $\ddot{Q}(\widehat{\btheta})$ is block-diagonal with blocks $\ddot{Q}_{11}(\widehat{\bbeta})$ and $\ddot{Q}_{22}(\widehat{\balpha})$, then we have, for any unit vector $\mathbf{d}$, \begin{eqnarray*}C_{f_{Q},\mathbf{d}}=C_{1,\mathbf{d}}(\hat{\bbeta})+C_{2,\mathbf{d}}(\hat{\balpha}),\end{eqnarray*} where $C_{1,\mathbf{d}}(\hat{\bbeta})=2\mathbf{d}^{\top}\Delta_{1\bomega_0}^{\top}(-\ddot{Q}_{11})^{-1}\Delta_{1\bomega_0}\mathbf{d}$ and $ C_{2,\mathbf{d}}(\hat{\balpha})=2\mathbf{d}^{\top}\Delta_{2\bomega_0}^{\top}(-\ddot{Q}_{22})^{-1}\Delta_{2\bomega_0}\mathbf{d},$ with
\begin{eqnarray*}\bDelta_{1\bomega_0}&=&\bDelta_{\bbeta}=\displaystyle\frac{\partial^2}{\partial\bbeta\partial\bomega^{\top}}\left\{Q(\btheta,\bomega|\widehat{\btheta})\right\}|_{\bomega_0}~~~\textrm{and} \\ \bDelta_{2\bomega_0}&=&\bDelta_{\balpha}=\left(\bDelta^{\top}_{\alpha 1},\bDelta^{\top}_{\alpha 2},\bDelta^{\top}_{\alpha 3}\right)^{\top},\end{eqnarray*} with $\bDelta_{\alpha_k}=\displaystyle\frac{\partial^2}{\partial\alpha_k\partial\bomega^{\top}}\left\{Q(\btheta,\bomega|\widehat{\btheta})\right\}|_{\bomega_0}\,\,{\rm for} \,k=1,2,3$ where $\ddot{Q}_{11}$ \ and \ $\ddot{Q}_{22}$ are defined in Equations (\ref{hessian1}) and (\ref{hessian2}).}
\noindent $\bullet$ \emph{Response perturbation:} We replace $V_{i}$ by $V_{i}(\omega)= V_{i}+\omega_{i}$ where ${\bomega}=\left(\omega_1,...,\omega_n\right)^{\top}$ is an $n \times 1$ vector with $\omega_i>0$. Recall that $V_i=V_{0i}$ if $C_i=0$ and $V_i=[V_{1i},V_{2i}]$ if $C_i=1$. Now, substituting $V_{i}(\omega)$ into Equation (\ref{CensL1}), we can write the perturbed model as:
\begin{eqnarray}
C_i =
\left\{\begin{array}{ccc}
1 & if & V_{1i}(\omega)\leq Z_{i}(\omega)\leq V_{2i}(\omega)\,, \\
0 & if & Z_{i}(\omega) = V_{0i}(\omega)\,.
\end{array}\right. \nonumber
\end{eqnarray}
where $\V_{1}(\bomega)=\V_1+\bomega,~\V_2(\bomega)=\V_2+\bomega$ and $\Z(\bomega)=\Z-\bomega$. Hence, the perturbed Q-function $Q(\btheta|\widehat{\btheta},\bomega)$ is as in Subsection 2.2, with $\widehat{\Z}$ and $\widehat{\Z\Z^{\top}}$ being replaced by $\widehat{\Z_{\omega}}=\widehat{\Z}-\bomega$ and $\widehat{\Z_{\omega}\Z_{\omega}^{\top}}=\widehat{\Z\Z^{\top}}-\widehat{\Z}\bomega^{\top}-\bomega\widehat{\Z}^{\top}+\bomega\bomega^{\top}$, respectively.
Under this perturbation scheme, the vector of no perturbation is given by $\bomega_0=\left(0,\ldots,0\right)^{\top}$ and $\bDelta_{\bomega_0}$ has the following elements:
\begin{equation}\label{P}
\bDelta_{\bbeta} = \X^{\top}\bSigma^{-1}\,\,\,\mbox{ and } \,\,\,
\bDelta_{\alpha_k}=-\widehat{\Z}^{\top}\left(\I_n-2\bP\right)^{\top}\frac{\partial\bSigma^{-1}}{\partial\alpha_k},\,\,\,k=1,2,3,
\end{equation}
where $\widehat{\Z}=(\widehat{Z}_1,\ldots,\widehat{Z}_n)^{\top}$, $ \bP=\X(\X^{\top}\bSigma^{-1}\X)^{-1}\X^{\top}\bSigma^{-1}$ and $\displaystyle\frac{\partial\bSigma^{-1}}{\partial\alpha_k}$ is as defined in Subsection \ref{HQmatrix}.
\noindent $\bullet$ \emph{Scale matrix perturbation:} In order to study the effects of perturbation of the scale matrix, we consider the perturbation scheme of the form $\bSigma(\bomega)=\bD(\bomega)\bSigma$, where $\bD(\bomega)$ is an $n \times n$ diagonal matrix with value $\omega_i$ in {\it i-}th diagonal element. Under this scheme, the non-perturbed model is obtained when $\bomega_0=\left(1,\ldots,1\right)^{\top}$. Thus, considering this perturbation scheme, $\bDelta_{\bomega_0}$ is a $(p+3)\times n $ matrix and has components given by:
\begin{equation*}\bDelta_{\bbeta i}=-\X^{\top}\bSigma^{-1}\bd(i)\left(\I_n-\bP\right)\Z~~~~\textrm{and}~~~\bDelta_{\balpha} = \PR{\Delta_{\alpha_{k i}}}\end{equation*} where $\Delta_{\alpha_{k i}} =
-\frac{1}{2}\chav{\textrm{tr}\PR{\widehat{\Z\Z^{\top}}\bSigma^{-1}\frac{\partial\bSigma}{\partial\alpha_k}\bSigma^{-1}\bd(i)}-\widehat{\Z}^{\top}\bSigma^{-1}\bP\frac{\partial\bSigma}{\partial\alpha_k}\bSigma^{-1}\bd(i)(2\I_n-\bP)\widehat{\Z}}, $
for $k=1,2,3$ and $\ii$. $\bd(i)$ is a $n \times n$ matrix with {\it i-}th diagonal element equal to one and the others equal to zero, $\displaystyle\frac{\partial\bSigma^{-1}}{\partial\alpha_k}$ and $\bP$ are as defined in Subsection \ref{HQmatrix} and Equation (\ref{P}), respectively. \\
\noindent $\bullet$ \emph{Explanatory variable perturbation:} We replace $\X$ by $\X(\bomega)=\X+\bW$ in the perturbed Q-function, with $\bW=\bomega\mathbf{1}^{\top}$ where $\bomega=\left(\omega_1,\dots,\omega_n\right)^{\top}$ and $\mathbf{1}$ is a $p\times 1$ vector of ones, so $\bW$ is an $n\times p$ matrix. Thus, considering the non-perturbed vector $\bomega_0=\left(0,\ldots,0\right)^{\top}$, the $\bDelta_{\bomega_0}$ has the following elements:
\begin{equation*}\bDelta_{\bbeta i}=\widehat{\Z}^{\top}(\I_n-2\bP)^{\top}\bSigma^{-1}\bW_i ~~~~\textrm{and}~~~\bDelta_{\balpha} = \PR{\Delta_{\alpha_{k i}}},\end{equation*} where ${\Delta_{\alpha_{k i}}} =
\widehat{\Z}^{\top}(\I_n-2\bP)^{\top}\frac{\partial\bSigma^{-1}}{\partial\alpha_k}\bW^{(1)}_i(\X^{\top}\bSigma^{-1}\X)^{-1}\X^{\top}\bSigma^{-1}\widehat{\Z},$
for $k=1,2,3$ and $\ii$. $\bW^{(1)}_i$ is a $n \times p$ matrix with the {\it i-}th row equal to one and the others equal to zero, $\displaystyle\frac{\partial\bSigma^{-1}}{\partial\alpha_k}$ and $\bP$ are as defined in Subsection \ref{HQmatrix} and Equation (\ref{P}), respectively.\\
{\bf Note also that it is not possible to give details for all the perturbation schemes that are of interest. However, as long as we can find an appropriate $\bomega$, and as long as the perturbed complete
data log-likelihood function $\ell_c(\btheta|\bC,\bV,\Z)$ is smooth enough, so that the required derivatives in the diagnostic measures are well defined, we can conduct the local influence analysis without much difficulty}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Real world data analysis}\label{realdata}
To exemplify the method developed in this work, we considered two dataset: \emph{depths of a geological horizon}, previously analyzed by \cite{dubrule1986interpolation} and \cite{de2005bayesian}, where the observations represent depths of a geological horizon measured at 100 locations and the \emph{Missouri dioxin contamination site}, analyzed by \cite{zirschky1986geostatistical}, measuring the level of contamination by dioxin at sampled points around the road. The computational procedures were implemented using the R software \citep{rmanual}, through the package {\it CensSpatial} \citep{CensSpatial.2016}.
\subsection{\textbf{Depths of a geological horizon}}
In this dataset, previously analyzed by \cite{dubrule1986interpolation} and \cite{de2005bayesian}, the observations are spread over a region of about 9 by 5 km and represent depths of a geological horizon measured at 100 locations, where 69 points are fully observed and 31 points are censored points, which are divided into left- and right- censored points. Figure \ref{dados} (a) shows a schematic description of the region displaying the sampling locations and the depth measurements. It is important to note that a "$\bullet$" represents an exact observed value, while "$\blacktriangledown$" represents a left-censored observation, and "$\blacktriangle$" indicates a location with a right-censored observation. The numbers above the symbols are the observed value for $C_i=0$ and the limit of detection for $C_i=1$, e.g., $\overset{1000}{\blacktriangle}$ means that at location $i$, $V_i=[1000,\infty)$ is observed.
\begin{figure}[ht]
\centering
\subfigure[]{\includegraphics[scale=0.5]{dadosDepth.eps}}
\subfigure[]{\includegraphics[scale=0.5]{variog.eps}}
\caption[]{\label{dados} Depth dataset. (a) Depths of a geological
horizon measured at each location. $\bullet$ represents an exact
observed value, $\blacktriangledown$ represents a left censored value and
$\blacktriangle$ represents a right censored value, the values above the
symbols indicate observed value for uncensored points and limit of
detection for censored points. (b) Plot of empirical variogram represented by $\circ$ and
theoretical variogram represented by solid line.}
\end{figure}
The depth data were transformed and their original units remain unknown for confidentiality reasons. For additional details about this dataset, we refer to \cite{de2005bayesian}.\\
\noindent \emph{(a) Model specification and preliminary analysis}\label{exploratory}
We consider the model:
$$Z_i=\mu+\epsilon_i,$$
with different correlation functions for the stochastic errors $\epsilon_i$, $i=1,\ldots, 100$. In this case, we have the presence of left and right censoring, simultaneously. Thus, in order to estimate the parameters of the model, the SAEM algorithm for censored data was implemented as described in Subsection \ref{SAEMspatial}. We choose a Monte Carlo sample size of $m=20$, a maximum number of iterations $W=150$ and a cutoff point $c=0.2$. {In our empirical studies, different choices of these arguments provided the same results for the SAEM estimates}. The computational procedures were implemented using the R software \citep{rmanual} and the codes are available from the authors upon request. To obtain information about the proper correlation function, as well as the initial values of the parameters that will be used in the SAEM algorithm, in Figure \ref{dados} (b) we depict the theoretical variogram using exponential, Gaussian and Mat\'{e}rn family with $\kappa=0.2$ and $\kappa=1$ correlation functions, along with the empirical variogram. From this figure, it can be noted that the first plotted ordinate is approximately 10, suggesting a relative nugget variance i.e., depth has some measurement error variance. The rising curve of sample variogram ordinates, leveling out at a distance of around 5, corresponds to a positive spatial correlation decaying with distance, and the sill is close to 30. This information suggests that the strength of spatial association is moderate. Additionally, it can be seen that the theoretical variogram of the exponential covariance function, is the closest to the empirical variogram.
\begin{table}[h!]
\begin{center}
\caption{Depth dataset. ML estimation under Exponential correlation function.}\label{ajusteDepth} \vskip 3mm
\begin{tabular*}{\textwidth}{l @{\extracolsep{\fill}} cccccccc}
\toprule
& $\widehat{\mu}$ & $\widehat{\sigma}^2$ & $\widehat{\rho}$ &
$\widehat{\tau}^2$ & &loglik & AIC & BIC \\
\midrule
Exponential & 1001.72 & 7.08 & 6.21 & 5.34 & & -165.87 & 339.74 & 350.16\\
\bottomrule
\end{tabular*}
\end{center}
\end{table}
As we are interested in presenting an diagnostic analysis, we proceed with estimation using the initial values described above and using the exponential correlation function. The ML estimates are presented on Table \ref{ajusteDepth}. Other details on the estimation and interpretation of the parameter estimates are omitted for brevity.\\
\noindent \emph{(b) Diagnostic analysis}
First, we consider the local influence analysis, based on $M(0)$, with interest focusing on $\btheta$ and using the following perturbation schemes: response perturbation, explanatory variables perturbation and scale matrix perturbation. It is important to stress that in local influence analysis, there are
no general rules so far for selecting the benchmark. Hence, we follow the criterion suggested by \cite{matos2013influence}, i.e.,
$M(0)_i>\overline{M}(0)+ 3SM(0)$, $i=1,\ldots,100,$ to discriminate whether an observation is influential or not.
Figure \ref{Localdepth} \subref{fig:subfig1}, \subref{fig:subfig2} and \subref{fig:subfig3} shows the index plot of $M(0)$ under the three perturbation schemes. We find that subjects $\#31,\#33,\#35$ and $\#94$ appear as influential under response perturbation and explanatory variable schemes. Since we do not have explanatory variables, the influential points showed by this scheme can be viewed as influential points in the expected average of the variable response $\mu$. For perturbation on $\bSigma$, we find that observations $\#14,\#34$ and $\#81$ appear as influential, which may indicate a more significant impact on the neighbors than the others points.
\begin{figure}[ht]
\centering
\subfigure[]{\includegraphics[scale=0.35,trim = 0mm 0mm 0mm 15mm,clip]{M0zdepth.eps}\label{fig:subfig1}}~
\subfigure[]{\includegraphics[scale=0.35,trim = 0mm 0mm 0mm 15mm,clip]{M0xdepth.eps}\label{fig:subfig2}}~
\subfigure[]{\includegraphics[scale=0.35,trim = 0mm 0mm 0mm 15mm,clip]{M0sigdepth.eps}\label{fig:subfig3}}\\
\subfigure[]{\includegraphics[scale=0.35,trim = 0mm 0mm 0mm 15mm,clip]{GDdepth2D.eps}\label{fig:subfig4}}~
\subfigure[]{\includegraphics[scale=0.35,trim = 0mm 0mm 0mm 15mm,clip]{QDdepth.eps}\label{fig:subfig5}}
\caption[]{\label{Localdepth}Depth dataset. Index plot of $M(0)$, using
exponential covariance function, under three perturbation schemes: \subref{fig:subfig1} response
perturbation; \subref{fig:subfig2} exploratory variable perturbation and
\subref{fig:subfig3} scale matrix perturbation. Index plot of \subref{fig:subfig4} $GD_i$ and
\subref{fig:subfig5} $QD_i$ using the exponential covariance function.}
\end{figure}
Finally, in order to evaluate the effect on the ML estimates when some observations are deleted, we analyze the $GD_i(\btheta)$ and $QD_i(\btheta)$ distance, which is depicted in Figure \ref{Localdepth} \subref{fig:subfig4} and \subref{fig:subfig5}. The plot reveals that case $\#34$ is potentially influential on the parameter estimates.\\
\noindent \emph{(c) The impact of the detected influential observations}
From the diagnostic analysis (global influence and local influence), we found that observations $\#14,\#31,~$ $\#33,\#34,\#35,\#81$ and $\#94$ are potentially influential. Thus, in order to reveal the impact of these observations, on the parameter estimates, we refitted the model individually eliminating each one of these seven cases. Table \ref{RC1} shows the relative changes (in percentage) of each parameter estimate defined by:
\begin{eqnarray}\label{RCApl1}
RC_{\widehat{\gamma}}
&=&\left|\frac{\widehat{\gamma}-\widehat{\gamma}_{[i]}}{\widehat{\gamma}}\right|,
\end{eqnarray}
where $\widehat{\gamma}=\widehat{\beta},\widehat{\sigma^2},\widehat{\rho}$ or $\widehat{\tau^2}$ and $\widehat{\gamma}_{[i]}$ denotes the ML estimate of $\widehat{\gamma}$ after the $i$th observation of ${\bf Z}$ is removed.
Note that the influential points identified by the scale matrix perturbation ($\#14,\#34$ and $\#81$) show significant changes in the parameter estimate, indicating the need for special attention in the estimating of $\Sigma$.
\begin{table}[h!]
\begin{center}
\caption{\label{RC1} Depth dataset. Relative changes [RC (in \%)].}
\small
\begin{tabular*}{\textwidth}{c @{\extracolsep{\fill}} ccccccc}
\toprule
Dropped & $\PR{\#14}$ & $\PR{\#31}$ & $\PR{\#33}$ & $\PR{\#34}$ & $\PR{\#35}$ & $\PR{\#81}$ & $\PR{\#94}$ \\
Coordinate & (5.44,1.54) & (8.18,0.18) & (0.27,0.49) & (0.43,1.59) & (0.65,4.26) & (6.44,3.71) & (8.03,4.36) \\
\midrule
$RC_{\hat{\mu}}$ & 0.0288 & -0.0163 & 0.0076 & 0.0042 & 0.0067 & 0.0051 & 0.0157 \\
$RC_{\hat{\sigma^2}}$ & 7.9784 & 2.2849 & -1.0510 & 10.9432 & 1.1388 & 9.3227 & 2.2724 \\
$RC_{\hat{\rho}}$ & 38.2575 & -11.6565 & -22.4489 & 36.5394 & 2.0544 & 17.7969 & -8.8715 \\
$RC_{\hat{\tau^2}}$ & 17.2380 & -0.4075 & -3.0620 & 16.8800 & -1.6016 & 8.8114 & 0.1547 \\
\midrule
$MSPE$ & 1.5672 & 1.5500 & 1.7891 & 1.3059 & 1.7030 & 1.5848 & 1.4656 \\
$RC_{MSPE}$ & 0.0007 & 0.0102 & 0.1424 & 0.1661 & 0.0874 & 0.0119 & 0.0642 \\
\bottomrule
\end{tabular*}
\end{center}
\end{table}
It is also important to verify the impact that the influential observations may cause on the prediction, once it is fundamental on the geostatistical analysis and this is one of the attractive features of the SAEM algorithm. Thus, in order to evaluate the impact of influential points on prediction, we proceed with an experiment where $5\%$ of the data were randomly removed and used later for prediction purposes and then we measured in some way the distance between the real value and the predicted value. To assess the prediction quality, we consider only the fully observed values to be predicted, so the true observed values can be compared with the predicted values obtained. For instance, when $5\%$ of data are taken, it means that $5$ fully observed values are separated, with $31$ censored observations and $64$ fully observed values remaining, which will be used for estimation via the SAEM algorithm.
As a measure of prediction quality, the mean square prediction error (MSPE) is used to calculate the distance between the real value and the predicted value, which is defined by
\citep[see,][]{fridley2006data}:
\begin{eqnarray}\label{MSPE}
MSPE &=& \frac{\sum_{i=1}^n (Z_i-\widehat{Z}_i)^2}{n},
\end{eqnarray}
where $Z_i$ is the observed value, $\widehat{Z}_i$ is the predicted value and $n$ is the number of samples to be predicted. Note that the best predicted model is determined by the lowest value of MSPE.
Then, to analyze the impact of influential points on the prediction, we calculated the MSPE considering all observations and individually eliminating each of these seven influential cases. With these measures, it is possible to calculate the relative changes of MSPE considering all observations regarding the other MSPE, defined by
\begin{eqnarray}\label{RCofMSPE}
RC_{MSPE}
&=&\left|\frac{MSPE-MSPE_{[i]}}{MSPE}\right|.
\end{eqnarray}
In Table \ref{RC1}, we show the results of $MSPE$ and $RC_{MSPE}$ eliminating each one of the influential cases. It is important to note that the MSPE obtained with all influential points was $1.5661$.
%{\bf
Note that dropping the influential points to proceed with prediction, in a general context, is the same as losing information. But in spatial data, the relation between one point and its neighbors is very important and this loss of information does not always negatively affect the prediction. An example of this fact is the two points presented as the biggest changes ($\#33,\#34$) with relation to the MSPE. The removal of point $\#34$ improves the prediction, which may have occurred because of the big difference between the observed value of point $\#34$ and its neighbors. While, the neighborhood of point $\#33$ shows similar observed values to this point, its removal produces a loss of information and worse prediction regarding to the prediction of the influential point $\#33$. Figure C.1, provided in the Supplementary Material, shows the location of each influential observation.
%}
\subsection{\textbf{Missouri dioxin contamination site}}
In this section, we considered a dataset reported in \cite{zirschky1986geostatistical} with 127 observations distributed in an area of $3600m\times 65m$ on the shoulders of a country road located in Missouri, USA. The observations correspond to a level of contamination by dioxin (2,~3,~7,~8-tetrachlorodibenzo-p-dioxin or TCDD) at sampled points along the road. The spatial directions are the X-direction (measured in $1/100$ feet), representing direction parallel to the road and the Y-direction (measured in feet), representing the direction perpendicular to or away from the road. The road is located at the Y coordinate of 30. The shoulder of the road was divided into long transects in the X direction, most 200 feet, in which eight samples were taken. The eight samples were aggregated together to give one measurement per transect. For illustration purposes, we treat the values reported as coming from one sampled location, with the X coordinate indicating the start of the transect, see Figure \ref{VeroVsK}(a). Forty-three percent of the observations {(55 sites)} were censored, falling below some LOD. The level of detections range from 0.10 to 0.79 mg/kg.
\begin{figure}[!h]
\centering
\subfigure[]{\includegraphics[scale=0.5,trim = 27mm 0mm 20mm 15mm,clip]{dados.eps}}
\subfigure[]{\includegraphics[scale=0.35,trim = 27mm 0mm 40mm 15mm,clip]{LogLik.eps}}
\caption[]{\label{VeroVsK} Missouri data. (a) Proportional TCDD
observed on each location. $\circ$ represents an observed value and
$\bullet$ represents a censored value. (b) The profile
log-likelihood versus $\kappa$ using the Mat\'{e}rn family.}
\end{figure}
\noindent (a) \emph{Model specification and preliminary analysis}
We propose to fit the model
$$\log\left\{Z_i\right\}=\mu+\epsilon_i,$$
with different covariance function for the stochastic errors $\epsilon_i$, $i=1,\ldots, 127$. Considering the left censoring in the response variable, we implemented the SAEM algorithm, as explained in Subsection \ref{SAEMspatial}, to estimate the parameters. We chose a Monte Carlo sample size of $m=20$, a maximum number of iterations $W=150$ and a cutoff point $c=0.2$. These computational procedures were implemented using the R software \citep{rmanual}. The results of the ML estimates, using the spherical, exponential and Mat\'{e}rn covariance functions are presented in Table \ref{ajusteAIS}. Notice that, although the spherical covariance function is widely used in classic geostatistics, this structure does not present a good fit, showing the worst value for the log-likelihood and information criteria.
\begin{table}[h!]
\begin{center}
\caption{\label{ajusteAIS} Missouri data. ML estimates under different covariance functions.}
\small
\begin{tabular*}{\textwidth}{c @{\extracolsep{\fill}} ccccccc}
\toprule
&&&& & \multicolumn{3}{c}{Criteria} \\
\cmidrule{6-8}
& $\widehat{\mu}$ & $\widehat{\sigma^2}$ & $\widehat{\rho}$ & $\widehat{\tau^2}$ & loglik. & AIC & BIC \\
\midrule
Spherical & -1.4137 & 3.7597 & 0.2000 & 3.4526 & -216.2290 & 440.4580 & 451.8348 \\
{\bf Exponential} & {\bf -2.0122} & {\bf 4.8016} & {\bf 14.0557} & {\bf 0.2445} & {\bf -143.8896} & {\bf 295.7793} & {\bf 307.1560} \\
Mat\'{e}rn ($\kappa=0.75$) & -1.8044 & 4.4306 & 8.4704 & 0.3751 & -144.0814 & 296.1627 & 307.5395 \\
Mat\'{e}rn ($\kappa=1.00$) & -1.7030 & 4.1682 & 6.0301 & 0.4276 & -144.5400 & 297.0800 & 308.4567 \\
\bottomrule
\end{tabular*}
\end{center}
\end{table}
The Gaussian and exponential covariance functions are special cases of the Mat\'{e}rn family of covariance functions. In Figure \ref{VeroVsK} (b), we show the values of the log-likelihood using the Mat\'{e}rn covariance after applying the SAEM algorithm fixing $\kappa$ at values $0.1, 0.2, 0.3, \cdots, 2.0$. Thus, we choose $\kappa = 0.5$, which maximizes the profile log-likelihood and corresponds to the exponential covariance function. In the following we proceed with diagnostics analysis using the exponential covariance function (or Mat\'{e}rn ($\kappa = 0.5$)). As we currently focus on exploring influence diagnostics, details of the estimation and interpretation of the parameter estimates are omitted for brevity. \\
\noindent \emph{(b) Diagnostic analysis}
As the first application, we consider the local influence analysis, based on $M(0)$, with interest focusing on $\btheta$ and considering the following perturbation schemes: response perturbation, explanatory variables perturbation and scale matrix perturbation. The criterion $M(0)_i>\overline{M}(0)+ 3SM(0)$, $i=1,\ldots,127,$ was used to discriminate whether an observation is influential or not.
Figure \ref{Local} \subref{subfig:mis:1}, \subref{subfig:mis:2} and \subref{subfig:mis:3} shows the index plot of $M(0)$, under the three perturbation schemes. We find that subjects $\#40$, $\#42$, $\#45$, $\#47$ and $\#48$ appear as influential under response and explanatory variable perturbations. Since we do not have explanatory variables, the influential points showed by this scheme can be viewed as influential points in the expected average of the variable $\mu$. For perturbation on $\bSigma$, we find that observations $\#40$ and $\#45$ appear as influential, which may indicate a more significant impact on neighbors than the others points.
\begin{figure}[ht]
\centering
\subfigure[]{\includegraphics[scale=0.45,trim = 27mm 0mm 40mm 15mm,clip]{M0y.eps}\label{subfig:mis:1}}~
\subfigure[]{\includegraphics[scale=0.45,trim = 27mm 0mm 40mm 15mm,clip]{M0x.eps}\label{subfig:mis:2}}~
\subfigure[]{\includegraphics[scale=0.45,trim = 27mm 0mm 40mm 15mm,clip]{M0sig.eps}\label{subfig:mis:3}}\\
\subfigure[]{\includegraphics[scale=0.27]{Rplot_GDiMissouri.eps}\label{subfig:mis:4}}~
\subfigure[]{\includegraphics[scale=0.27]{Rplot_GDibetaMissouri.eps}\label{subfig:mis:5}}~
\subfigure[]{\includegraphics[scale=0.27]{Rplot_GDialphaMissouri.eps}\label{subfig:mis:6}}~
\caption[]{\label{Local}Missouri data. Index plot of $M(0)$ using exponential covariance and considering: \subref{subfig:mis:1} Response perturbation; \subref{subfig:mis:2} exploratory variable perturbation and \subref{subfig:mis:3} scale matrix perturbation. Index plot of \subref{subfig:mis:4} $GD_i$, \subref{subfig:mis:5} $GD({\bbeta})$ and \subref{subfig:mis:6} $GD({\balpha})$ for the TCDD data using exponential covariance function.}
\end{figure}
It is important to stress that one of the points indicated as influential, observation ($\#40$), is the maximum value of dioxin observed and it is located on the boarder of the road, where the dioxin was dumped, so it is expected to be an influential point. The other points are located in a perpendicular line (see Figure C.2 in the Supplementary Material) to the road which contains observation $\#40$. Considering the response and explanatory perturbations, which compare the observed data in relation its neighborhood, this may represent the direction where the dioxin was spread.
In order to evaluate the effect on the ML estimates when some observations are deleted, we analyze the $GD_i(\btheta)$ distance, which is depicted in Figure \ref{Local} \subref{subfig:mis:4}. The plot reveals that once again cases $\#40, \#45, \#47$ and $\#48$ are potentially influential on the parameter estimates. Figure \ref{Local} \subref{subfig:mis:5} and \subref{subfig:mis:6} present the index plot of $GD_i(\bbeta)$ and $GD_i(\balpha)$, respectively. From these figures, we see that the observations $\#40$, $\#45$, $\#47$ and $\#48$ are influential with regard to the parameter $\balpha$.\\
\noindent \emph{(c) The impact of the detected influential observations}
The diagnostic analysis indicated five observations $(\#40, \#42, \#45, \#47$ and $\#48)$ as potentially influential. In order to reveal the impact of these five observations, on the parameter estimates, we refitted the model individually eliminating each one of these five cases. Thus, in Table \ref{RC} we show the relative changes $(\%)$, of each parameter estimate, as presented in Equation (\ref{RCApl1}).
Note from Table \ref{RC} that significant changes are observed only for the nugget effect $\tau^2$, indicating the needed for special attention on estimation of the variance.
\begin{table}[htb]
\small
\caption{\label{RC} Missouri data. Relative changes, MSPE and relative changes of the MSPE [RC (in \%)].}
\begin{tabular*}{\textwidth}{c @{\extracolsep{\fill}} ccccc}
\toprule
Dropped & $\PR{\#40}$ & $\PR{\#42}$ & $\PR{\#45}$ & $\PR{\#47}$ & $\PR{\#48}$ \\
\midrule
$RC_{\hat{\mu}}$ & 1.7675 & 1.1246 & 0.1962 & 0.3180 & 1.0575 \\
$RC_{\hat{\sigma^2}}$ & 1.6641 & 0.5709 & 2.0432 & 0.7147 & 2.3936 \\
$RC_{\hat{\rho}}$ & 2.8503 & 4.3562 & 5.0879 & 5.342 & 2.3200 \\
$RC_{\hat{\tau^2}}$ & 18.7987 & 15.9953 & 22.5714 & 24.773 & 22.5531 \\
\midrule
$MSPE$ & 0.666 & 0.8871 & 0.8466 & 0.8436 & 0.8606\\
$RC_{MSPE}$ & 0.2146 & 0.0459 & 0.0018 & 0.0053 & 0.0146\\
\bottomrule
\end{tabular*}
\end{table}
To check the impact of influential points on prediction, we proceeded with the same experiment previously define, where $5\%$ of the data were randomly set aside to use in the prediction stage and then we used Equations (\ref{MSPE}) and (\ref{RCofMSPE}) to measure the quality of prediction and the impact of influential points, respectively. Table \ref{RC}, also shows this results, note that the MSPE obtained with all influential points was $0.8481$.
Note that the observation $\#40$ shows the biggest relative change but a better MSPE, what can be explained by the fact that this observation is the maximum value of dioxin observed and an expected influential value showing the attention needed for influential points on prediction.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{\color{blue}
\section{Simulation Studies}\label{simulation}
In order to assess the performance of the proposed diagnostic measures for identifying outliers, we consider a left censored spatial linear model (SCL) specified as:
\begin{equation}\label{sim:eq}
\Z=\beta_0+\beta_1\x_1+ \beta_2\x_2+\bepsilon,
\end{equation}
where $\bepsilon\sim N_n(\mathbf{0},\bSigma)$ with $\bSigma=[C(\sa_i,\sa_j)]=\tau^2\mathbf{I}_n+\sigma^2\R(\rho)$, considering different censoring levels and Mat\'{e}rn($\kappa=0.75$)) covariance function. We set $\bbeta^{\top}=(\beta_0,\beta_1,\beta_2)=(1,3,-1)$, $\sigma^2=3$, $\rho=3$, $\tau^2=2$, $n=100$ and $\x_i^{\top}=(1, x_{i1},x_{i2})$, where $x_{i1}\sim U(0,1)$ and $x_{i2}\sim U(2,5)$, for $i=1,\ldots,100$.
\subsection{\textbf{First simulation study}}
This study explores the accuracy of the proposed diagnostics measures in detecting two outliers under different levels of censoring ($10\%,20\%,30\%$), considering a Mat\'{e}rn covariance function with $\kappa=0.75$. In order to do this, we consider two influential points under two different scenarios: (1) The points are spatially-near to each other, and (2) the points are spatially-far. The points chosen for this study are $\#34$ and $\#44$, located at coordinates $(4,4),(4,5)$ for scenario 1, and $\#34$ and $\#78$ for scenario 2, located at the coordinates $(4,4),(8,8)$.
In this simulation we generated only one sample under three levels of censoring and the simulated samples were perturbed by replacing the selected observations by $z = z+b*sd(\zp)$, where $b$ assumes the values $2,4$ and $6$ . The larger is the value of $b$ the larger is the perturbed point. Thus, $b=2$ indicates low perturbation, $b=4$ indicates moderate perturbation and $b=6$ indicates high perturbation. For the sake of simplicity, only the results of $30\%$ of censoring are described, since similar results are obtained for other censoring levels. \\
\noindent {\it Scenario 1: Points spatially-near}
Figure C.3, provided in the Supplementary Material, depicted the simulated values of $z$ for each value of $b$. Note that for a low perturbation ($b=2$), the perturbed points are not far from the other points.
Following the approach described on Section 4, Figure \ref{FS1local} depicts the index plots of $M(0)$ for the response perturbation, explanatory variable perturbation and scale perturbation, respectively, along with the \cite{LeeXu04} benchmark computed for $c^{*}=3$. From this Figure we observe that for all values of $b$, at least one of the perturbed points is identified as an influential point. Note that under this scenario, each perturbed point has eight neighbors, with four neighbors in common and each other as a neighbor. Under this consideration, the local influence measures were less sensitive for low and moderate perturbation, which can be explained by the fact that the impact generated by one of the perturbed points was absorbed by the impact generated by the other one, i.e., the influence of both perturbed points over its neighbors can be similar, making harder to identify simultaneously the influential points.
\begin{figure}[!h]
\center \subfigure[ref1][]{\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0z30_2.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0x30_2.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0sig30_2.eps}}
\qquad \subfigure[ref2][]{\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0z30_4.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0x30_4.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0sig30_4.eps}}
\qquad \subfigure[ref3][]{\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0z30_6.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0x30_6.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0sig30_6.eps}}
\caption[]{\label{FS1local} Simulation Study 1. Scenario 1. Index plot of $M(0)$ for
(right) response perturbation; (middle) explanatory variable perturbation and
(left) scale matrix perturbation, using Mat\'{e}rn Covariance function with
$30\%$ of censoring and $b = 2(a),4(b),6(c)$.}
\end{figure}
In Figure C.4, provided in the Supplementary Material C, we present the index plot of the global influence measures, $GD_i$ and $QD_i$. Here we note similar conclusions than those under the local influence approach.\\ %However, for these %global influence measures the impact of the influential points over the parameters estimates, then these %results may imply that in the case of two perturbed points spatially close, one generates more impact over the %estimates than the other.\\
\noindent {\it Scenario 2: Points spatially-far}
Like the first scenario, Figure C.5 (Supplementary Material) presents the perturbed values $z$ for each value of $b$. Note that for a low perturbation ($b=2$), the perturbed points are not far from the other points.
Following the approach described in Section 4, Figure \ref{FS2local} depicts the index plots of $M(0)$ for the response perturbation, explanatory variable perturbation and scale perturbation, respectively. Note that with perturbed points being spatially far from each other, the local diagnostics measures may easily identify those points as influential. Regardless the case of low perturbation ($b=2$), where the perturbed points do not stand out of the others.
\begin{figure}[!h]
\center \subfigure[ref1][]{\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0z2_302.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0x2_302.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0sig2_302.eps}}
\qquad \subfigure[ref2][]{\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0z2_304.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0x2_304.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0sig2_304.eps}}
\qquad \subfigure[ref2][]{\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0z2_306.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0x2_306.eps}~
\includegraphics[scale=0.37,trim = 0mm 10mm 10mm 10mm,clip]{M0sig2_306.eps}}
\caption{\label{FS2local}Simulation Study 1. Scenario 2. Index plot of $M(0)$ for
(right) response perturbation; (middle) explanatory variable perturbation and
(left) scale matrix perturbation, using Mat\'{e}rn Covariance function with
$30\%$ of censoring and $b = 2(a),4(b),6(c)$.}
\end{figure}
In Figure C.6 (Supplementary Material) we present the index plot of the global influence measures, $GD_i$ and $QD_i$. As expected, once again the perturbed points were detected as influential.
\subsection{\textbf{Second simulation study}}
The second study is a Monte Carlo experiment, which shows the capacity of the method to detect atypical points. Here, we consider as in the first study two atypical points spatially-near (observations \#34, \#44) and two atypical points spatially-far (observations \#34, \#78). We also consider the direction of influence, i.e., both points are perturbed in order to assume large values. In this study, the levels of perturbation are $b=2,4,6$. Since we have a benchmark to decide which point is influential or not, the diagnostic measure $M(0)$ was computed for 500 simulated datasets under $30\%$ of censoring proportions and using Mat\'{e}rn (with $\kappa=0.75$) covariance function.
\begin{table}[!h]
\begin{center}
\caption{Simulation study 2. The values in the table denotes the \%
of correctly identifying the influential observations using response
perturbation ($\z$), explanatory variable perturbation ($\X$) and
scale matrix perturbation ($\bSigma$) from 500 simulated datasets
under the SCL model with Mat\'{e}rn correlation
structure and different levels of perturbation denoted by $b$ for Scenario 1 and Scenario 2.}\label{SS1} \vskip 3mm
\small
\begin{tabular*}{\textwidth}{ll @{\extracolsep{\fill}} | lccc | lccc }
\toprule
& & \multicolumn{4}{c}{Scenario 1} & \multicolumn{4}{c}{Scenario 2}\\
\cmidrule(r){3-6}
\cmidrule(r){7-10}
& & &\multicolumn{3}{c}{Scheme of perturbation} & & \multicolumn{3}{c}{Scheme of perturbation}\\
& & & $\Z$ & $\X$ & $\bSigma$ & & $\Z$ & $\X$ & $\bSigma$ \\
\midrule
\multirow{3}{*}{\textbf{$b=2$}} & & $\#34$ & 53.0 & 58.0 & 69.4 & $\#34$ & 67.8 & 65.6 & 73.6 \\
& & $\#44$ & 60.6 & 63.8 & 71.6 & $\#78$ & 70.8 & 65.8 & 74.8 \\
& & $\#34,\#44$ & 26.0 & 31.0 & 42.8 & $\#34,\#78$ & 42.2 & 37.2 & 48.8 \\
\midrule
\multirow{3}{*}{\textbf{$b=4$}} & & $\#34$ & 78.0 & 74.4 & 79.8 & $\#34$ & 90.2 & 84.8 & 82.8 \\
& & $\#44$ & 94.8 & 89.0 & 93.8 & $\#78$ & 98.4 & 93.4 & 94.0 \\
& & $\#34,\#44$ & 73.0 & 63.6 & 73.6 & $\#34,\#78$ & 88.6 & 78.2 & 76.8 \\
\midrule
\multirow{3}{*}{\textbf{$b=6$}} & & $\#34$ & 81.6 & 75.2 & 72.8 & $\#34$ & 94.2 & 89.4 & 79.0 \\
& & $\#44$ & 99.8 & 98.0 & 99.8 & $\#78$ & 100.0 & 99.6 & 100.0 \\
& & $\#34,\#44$ & 81.6 & 73.2 & 72.6 & $\#34,\#78$ & 94.2 & 89.0 & 79.0 \\
\bottomrule
\end{tabular*}
\end{center}
\end{table}
As mentioned above, we consider the following perturbation $z_{i}=\parallel z_{i}\parallel+b*sd(\zp)$, for $i=34,44$ and $78$. Table \ref{SS1} reports the percentage of times the measures correctly identifies $z$ as the most influential. As expected, the percentage of correctly detecting atypical observations increases for increasing level of perturbation $b$.
However, it can be noted that the capacity of the method to detect both influential observations jointly under low perturbation ($b=2$) is poor when compared with individual detections, which is in general have good performance.
%}
{\bf To examine and investigate the asymptotic properties of the SAEM estimates, a simulation study was conducted and the results are presented in Appendix D provided in the Supplementary Material. We also provided a simulation study to measure the impact of misspecification of the covariance structure (Appendix E).}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}\label{concl}
This article proposes influence diagnostic tools for detecting influential observations in the context of spatial censored linear models. It extends the recently published works by \cite{assumpccao2014analysis} and \cite{de2014influence} which consider estimation and diagnostics of spatial linear models. Our proposed method relies on the $Q$-function, the conditional expectation of the logarithm of the complete-data likelihood, which facilitates the theoretical development of the stochastic approximation of the EM algorithm (SAEM) to obtain the maximum likelihood estimates of model parameters and the development of influence diagnostic measures. Explicit expressions are obtained for the Hessian matrix $\ddot{\bQ}$ and for the matrix $\Delta$ under different perturbation schemes. A simulation study compares the outlier detection accuracy under different censoring and perturbation schemes. For practical demonstration, the method is applied to a data from dioxin contaminated sites in Missouri and the proposed methods are implemented using the R software (codes available upon request from the first author), providing practitioners with a convenient tool for further applications in their domain.
{\bf As pointed out by an anonymous referee, there is an important issue that could not be handled by the \cite{cook86} or \cite{zhu2001case} approaches. It is related to the appropriate choice of a perturbation vector. As arbitrarily perturbing a model may lead to misleading inference about the cause of an influential effect, this issue is important in local influence analysis. Recently, \cite{zhu2007perturbation} proposed to use a metric tensor to select an appropriate
perturbation to a statistical model, extending the work of \cite{cook86} and \cite{zhu2001case}. However, there are many issues related to this method, such as the calculation of the influence measures and metric tensor
under different situations. An extension of this method to the SCL model is of great interest but
is beyond the scope of this paper.}
Future extensions of the work include the use of scale mixtures of normal distributions to accommodate heavy-tailed features, or the development of some diagnostics and tests for the model. Bayesian influence diagnostics, in the context of spatial censored linear models, can be treated via the Kullback--Leibler divergence, as proposed by \cite{Cancho2010}. Other extensions of the current work include, for example, diagnostics analysis in censored spatial data with measurement errors \citep{li2009spatial}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgements}
The authors are grateful to the Editor, Associate Editor and two referees for their helpful comments on an earlier version of this paper. The research of V\'{\i}ctor H. Lachos was supported by Grant 306334/2015-1 from Conselho Nacional de Desenvolvimento Cient\'{\i}fico e Tecnol\'{o}gico (CNPq-Brazil) and by Grant 2014/02938-9 from Funda\c{c}\~{a}o de Amparo \`{a} Pesquisa do Estado de S\~{a}o Paulo (FAPESP-Brazil). The research of Thais S. Barbosa was supported by CAPES. Larissa A. Matos acknowledges support from FAPESP-Brazil (Grant 2016/05420-6).
The research of Aldo M. Garay was supported by Grant 420082/2016-6 from CNPq-Brazil. This paper was written while V\'{\i}ctor H. Lachos was a visiting professor in the Department of Statistics at the University of Connecticut, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{wb_env}
\bibliography{biblio}
\label{lastpage}
\end{document}