APCTP Pre2022-0xx

#  A model explaining the new CDF II W boson mass linking to muon $g-2$ and
dark matter

Keiko I. Nagao<EMAIL_ADDRESS>Okayama University of Science, Faculty of
Science, Department of Applied Physics, Ridaicho 1-1, Okayama, 700-0005, Japan
Takaaki Nomura<EMAIL_ADDRESS>College of Physics, Sichuan University,
Chengdu 610065, China  Hiroshi Okada<EMAIL_ADDRESS>Asia Pacific
Center for Theoretical Physics, Pohang 37673, Republic of Korea  Department of
Physics, Pohang University of Science and Technology, Pohang 37673, Republic
of Korea

###### Abstract

We propose a model to explain the W boson mass anomaly reported by CDFII
collaboration that would suggest new physics. We introduce exotic fermions;
one isospin doublet vector-like lepton, one isospin singlet singly-charged
vector-like lepton and an isospin doublet inert scalar. These new fermions
also provide us with sizable anomalous magnetic moment (muon $g-2$) due to no
chiral suppression and nonzero mass difference between the real and imaginary
parts of neutral inert scalar bosons. This mass difference and doublet vector-
like leptons affect the oblique parameters, especially $\Delta T$, whose shift
from zero are the direct consequence of this W boson mass anomaly. We search
for the allowed parameter region to explain both the phenomenologies. In
addition, we discuss a dark matter candidate assuming the real part of the
inert scalar boson to be the one. We find that lighter DM mass is favored.

##  I Introduction

Recently, CDFII collaboration reported an updated result of the standard model
charged-gauge boson (W boson) mass CDF:2022hxs

$\displaystyle m_{W}=(80.433\pm 0.0064_{\rm stat}\pm 0.0069_{\rm syst})\ {\rm
GeV},$ (1)

which deviates from the standard model (SM) prediction by 7$\sigma$, where the
SM tells us $m_{W}=(80.357\pm 0.006)\ {\rm GeV}$. This disagreement also
appeared in the previous global combination of data from LEP, CDF, D0, and
ATLAS where they give the mass range of $m_{W}=(80.379\pm 0.012)\ {\rm GeV}$
ParticleDataGroup:2018ovx . This anomaly suggests new physics (NP) beyond the
SM Fan:2022dck ; Lu:2022bgw ; Athron:2022qpo ; Yuan:2022cpw ; Strumia:2022qkt
; Yang:2022gvz ; deBlas:2022hdk ; diCortona:2022zjy ; Du:2022pbp ;
Tang:2022pxh ; Blennow:2022yfm ; Cacciapaglia:2022xih ; Sakurai:2022hwh ;
Fan:2022yly ; Liu:2022jdq ; Lee:2022nqz ; Cheng:2022jyi ; Song:2022xts ;
Bagnaschi:2022whn ; Paul:2022dds ; Bahl:2022xzi ; Asadi:2022xiy ;
DiLuzio:2022xns ; Athron:2022isz ; Gu:2022htv ; Heckman:2022the ; Babu:2022pdn
; Heo:2022dey ; Du:2022brr ; Cheung:2022zsb ; Crivellin:2022fdf ; Endo:2022kiw
; Biekotter:2022abc ; Balkin:2022glu ; Krasnikov:2022xsi ; Ahn:2022xeq ;
Han:2022juu ; Zheng:2022irz ; Kawamura:2022uft ; Ghoshal:2022vzo ;
Perez:2022uil , and can be interpreted as the deviation of oblique parameters
Peskin:1991sw , especially $\Delta T$. In fact, the oblique parameters are
zero in the SM.

In this letter, we explain this anomaly by introducing exotic fermions; one
isospin doublet vector-like fermion $(L^{\prime})$, one isospin singlet
singly-charged vector-like fermion ($E$) and an isospin doublet inert scalar
($\eta$). These new fermions provide us with sizable muon $g-2$ due to no
chiral suppression. On the other hand, the doublet fermion negatively
contributes to $\Delta T$, while the doublet boson can positively contributes
to it. The current deviation tends to shift to the positive side of $\Delta
T$, and the doublet scalar is definitely required. We also show that these new
fields simultaneously connect to muon anomalous magnetic moment (muon $g-2$ or
$\Delta a_{\mu}$) and dark matter (DM) both of which would imply NP.

This letter is organized as follows. In Sec. II, we review our model and show
masses and mixings for new fields. In Sec. III, we formulate our oblique
parameters, muon $g-2$, and DM. In Sec.IV, we carry out the numerical analysis
and demonstrate our allowed region. Sec. V is devoted to our summary and
conclusion.

##  II Model setup and Constraints

|  $L_{L}$ |  $\ell_{R}$ |  $L^{\prime}$ |  $E$ |  $H$ |  $\eta$
---|---|---|---|---|---|---
$SU(3)_{C}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$SU(2)_{L}$ | $\bm{2}$ | $\bm{1}$ | $\bm{2}$ | $\bm{1}$ | $\bm{2}$ | $\bm{2}$
$U(1)_{Y}$ | $-\frac{1}{2}$ | $-1$ | $-\frac{1}{2}$ | $-1$ | $\frac{1}{2}$ | $\frac{1}{2}$
$Z_{2}$ | $+$ | $+$ | $-$ | $-$ | $+$ | ${-}$
Table 1:  Charge assignments of the new fields
$L^{\prime}\equiv[N^{\prime},E^{\prime}]^{T}$, $E$, and $\eta$ under
$SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times Z_{2}$. Here $L^{\prime}$ and
$E$ are vector-like fermions. $\eta$ is an inert scalar boson,and its real
neutral component is considered as a DM candidate.

In this section, we define our model. We introduce a vector-like exotic lepton
$L^{\prime}\equiv[N^{\prime},E^{\prime}]^{T}$ (weak isospin doublet) and a
singly charged-lepton $E$ (weak isospin singlet), and an inert doublet scalar
$\eta$ whose neutral component can be a DM candidate. Furthermore, we impose
odd parity under $Z_{2}$ discrete symmetry for new fields that assure the
stability of DM. The SM Higgs is denoted by $H$ and its VEV is given by
$[0,v/\sqrt{2}]^{T}$. Field contents and their charge assignments are
summarized in Table 1. Then, the renormalizable valid Lagrangian is given by

$\displaystyle-\mathcal{L}=$ $\displaystyle
y_{\ell_{ij}}\overline{L_{L_{i}}}H\ell_{R_{j}}+y_{L^{\prime}_{i}}\overline{E_{R}}\eta^{\dagger}L_{L_{i}}+y_{S_{i}}\overline{\ell_{R_{i}}}\eta^{\dagger}L^{\prime}_{L}+y_{E}\overline{L^{\prime}_{L}}HE_{R}+M_{E}\overline{E}E+M_{L^{\prime}}\overline{L^{\prime}}L^{\prime}+{\rm
h.c.},$ (2)

where $i,j=e,\mu,\tau$ are the flavor indices and we have abbreviated kinetic
terms for simplicity.

The scalar potential in our model is given by

$\displaystyle V$
$\displaystyle=-M_{H}^{2}H^{\dagger}H+M^{2}_{\eta}|\eta|^{2}+\lambda_{H}(H^{\dagger}H)^{2}+\lambda_{\eta}|\eta|^{4}+\lambda_{H\eta}(H^{\dagger}H)(\eta^{\dagger}\eta)$
$\displaystyle+\lambda^{\prime}_{H\eta}(H^{\dagger}\eta)(\eta^{\dagger}H)+\lambda^{\prime\prime}_{H\eta}[(H^{\dagger}\eta)^{2}+{\rm
h.c.}].$ (3)

We write $H$ and $\eta$ to be

$\displaystyle H=\begin{pmatrix}w^{+}\\\
\frac{1}{\sqrt{2}}(v+{h}+iG_{Z})\end{pmatrix},\
\eta=\begin{pmatrix}\eta^{+}\\\
\frac{1}{\sqrt{2}}({\eta_{R}}+i\eta_{I})\end{pmatrix},$ (4)

where ${h}$ corresponds to CP-even physical scalar boson state while $w^{+}$
and $G_{Z}$ are the Nambu-Goldstone(NG) bosons which are absorbed by the SM
weak gauge bosons $Z,~{}W$. Since we consider the neutral component of $\eta$
to be DM $\eta$ has zero VEV and no mixing with $H$; an inert scalar boson.
Here, we simply define these masses to be $m_{R},\ m_{I},\ m_{\eta^{\pm}}$,
respectively. The mass differences among inert scalar boson masses are arisen
from terms with couplings $\lambda^{\prime}_{H\eta}$ and
$\lambda^{\prime\prime}_{H\eta}$ such that

$m_{R}^{2}-m_{I}^{2}=\lambda^{\prime\prime}_{H\eta}v^{2},\quad
m^{2}_{\eta^{\pm}}-m^{2}_{R}=-(\lambda^{\prime}_{H\eta}+\lambda^{\prime\prime}_{H\eta})v^{2}.$
(5)

The difference between $m_{R}$ and $m_{I}$ is important to DM and muon $g-2$
as discussed later.

Exotic singly-charged fermions $(E,E^{\prime})$ mix each other after
spontaneous electroweak symmetry breaking, and the corresponding mass matrix
is written by

$\displaystyle{\cal M}_{E}=\left(\begin{array}[]{cc}M_{E}&m_{E}\\\
m_{E}&M_{L^{\prime}}\\\ \end{array}\right),$ (8)

where $m_{E}\equiv y_{E}v/\sqrt{2}$. The mass matrix can be diagonalized by a
unitary matrix $V_{E}$ as ${\rm diag}(M_{1},M_{2})=V_{E}^{\dagger}{\cal
M}_{E}V_{E}$, and $V_{E}$ is given by

$\displaystyle V_{E}=\left(\begin{array}[]{cc}c_{c}&-s_{c}\\\ s_{c}&c_{c}\\\
\end{array}\right),\quad\tan 2\theta_{c}=\frac{2m_{E}}{M_{E}-M_{L^{\prime}}},$
(11)

where $c_{c}(s_{c})$ is the shorthand notation of
$\cos\theta_{c}(\sin\theta_{c})$. We write mass eigenstates as $E_{1}$ and
$E_{2}$ whose masses are denoted as $M_{1}$ and $M_{2}$.

##  III Phenomenological Constraints and implications

In this section we summarize formulas for oblique parameters, muon $g-2$, SM
Higgs diphoton decay ratio and DM annihilation cross section via Yukawa
interactions.

###  III.1 Oblique parameters

The vector-like fermions $L^{\prime}(E)$ and an inert scalar doublet $\eta$
contribute to oblique parameters via vacuum polarization diagrams for
electroweak gauge bosons. Here, we compute $S$\- and $T$-parameters that are
influential in the W boson anomaly. They are given as follows: Peskin:1990zt ;
Ko:2021lpx

$\displaystyle\alpha_{em}S=4e^{2}\left(\sum_{\rm
VLF}\left[\frac{d}{dq^{2}}\Pi_{33}^{\rm VLF}-\frac{d}{dq^{2}}\Pi_{3Q}^{\rm
VLF}\right]_{q^{2}=0}\right),$ (12)
$\displaystyle\alpha_{em}T=\frac{e^{2}}{s_{W}^{2}c_{W}^{2}m_{Z}^{2}}\left(\sum_{\rm
VLF}\left[\Pi^{\rm VLF}_{\pm}(q^{2})-\Pi^{\rm
VLF}_{33}(q^{2})\right]_{q^{2}=0}\right),$ (13)

where the superscript VLF represents contribution from different combination
of our vector-like fermions $\\{E_{1},E_{2},N^{\prime}\\}$, and $\Pi^{{\rm
VLF}}(33,3Q,\pm)$ are contributions of beyond the SM. Also, $\Pi_{33,3Q,\pm}$
can be obtained from vacuum polarization contributions listed in Appendix A.
The contributions of $\eta$ to $S,\ T$ are simply given by Barbieri:2006dq

$\displaystyle\Delta S$
$\displaystyle\approx\frac{1}{2\pi}\int_{0}^{1}dxx(1-x)\ln\left[\frac{xm_{R}^{2}+(1-x)m_{I}^{2}}{m_{\eta^{\pm}}^{2}}\right],$
(14) $\displaystyle\Delta T$
$\displaystyle\approx\frac{1}{24\pi^{2}\alpha_{\rm
em}v^{2}}(m_{\eta^{\pm}}-m_{I})(m_{\eta^{\pm}}-m_{R}).$ (15)

In the numerical analysis, the deviation from the SM is discussed.

###  III.2  Muon anomalous magnetic dipole moment

New results on the muon $(g-2)$ were recently published by the E989
collaboration at Fermilab Muong-2:2021ojo :

$\displaystyle a^{\rm FNAL}_{\mu}=116592040(54)\times 10^{-11}.$ (16)

Combined with the previous BNL result, this means that the muon $(g-2)$
deviates from the SM prediction by 4.2$\sigma$ level Muong-2:2021ojo :

$\displaystyle\Delta a^{\rm new}_{\mu}=(25.1\pm 5.9)\times 10^{-10},$ (17)

and it could be a telltale sign of the physics beyond the SM.

Figure 1: An one loop diagram dominantly contributing muon $g-2$

In our model, the new dominant contribution to muon $g-2$ comes from the one
loop diagrams involving the vector-like leptons and neutral component of
$\eta$. The dominant contribution is given by the diagram in Fig. 1 and
consists of a product of $y_{L^{\prime}_{\mu}}$ and $y_{S_{\mu}}$ 111The
contributions proportional to $|y_{L^{\prime}}|^{2}$ or $|y_{S}|^{2}$ cannot
be dominant because of the chiral suppression. Thus, we neglect those terms
for simplicity.. Calculating the diagram muon $g-2$ is given by

$\displaystyle\Delta a^{\rm new}_{\mu}$
$\displaystyle=\frac{m_{\mu}}{(4\pi)^{2}}s_{c}c_{c}y_{L^{\prime}_{\mu}}y_{S_{\mu}}\sum_{a=1,2}\frac{(-1)^{a-1}}{M_{a}}\left[{F(q_{a},r_{Ra})}-{F(q_{a},r_{Ia})}\right],$
(18) $\displaystyle F(q,r)$
$\displaystyle\equiv\int_{0}^{1}dx\int_{0}^{1-x}dy\frac{1-2y}{1-x+q^{2}(x^{2}-x)+r^{2}x}\approx\frac{-1+r^{4}-2r^{2}\ln(r^{2})}{2(-1+r^{2})^{3}}+{\cal
O}(q^{2}),$ (19)

where $q_{i}\equiv\frac{m_{\mu}}{M_{i}}(<<1)$ and
$r_{ai}\equiv\frac{m_{a}}{M_{i}}$ ($a=R,\ I$, $i=1,2$). As can be seen from
the above formula, our new contribution exactly vanishes when $m_{R}=m_{I}$.
Thus, the (not so small) mass difference is crucially important to obtain the
sizable muon $g-2$. In our numerical analysis, we take the experimental region
at 3$\sigma$ interval.

###  III.3 Modification of $h\to\gamma\gamma$ decay width

New charged particles in the model modify the diphoton decay width of the SM
Higgs boson. Relevant interactions are given by

$\mathcal{L}\supset\frac{y_{E}}{\sqrt{2}}h(\overline{E^{\prime}_{L}}E_{R}+\overline{E}_{R}E^{\prime}_{L})+\lambda_{H\eta}vh\eta^{+}\eta^{-},$
(20)

where $E$ and $E^{\prime}$ are written by mass eigenstates through mixing Eq.
(11). We obtain deviation from the SM prediction from one loop diagrams in
which $\\{\eta^{\pm},E_{1},E_{2}\\}$ propagate inside the loop, such that

$\displaystyle\mu_{\gamma\gamma}\equiv\frac{\Gamma(h\to\gamma\gamma)}{\Gamma(h\to\gamma\gamma)_{\rm
SM}}\simeq\left|1+\frac{C_{E}\left(\frac{A_{1/2}(\tau_{1})}{M_{1}}-\frac{A_{1/2}(\tau_{2})}{M_{2}}\right)+\frac{\lambda_{H\eta}v^{2}}{2m_{\eta^{\pm}}^{2}}A_{0}(x_{\eta^{\pm}})}{A_{1}(\tau_{W})+\frac{4}{3}A_{1/2}(\tau_{t})}\right|^{2},$
(21) $\displaystyle
C_{E}=\frac{s^{2}_{c}c^{2}_{c}(M_{E}-M_{L^{\prime}})}{c^{2}_{c}-s^{2}_{c}},\
A_{1/2}(x)=-2[x+(1-x)f(x)],\ A_{1}(x)=2+3x+3(2x-x^{2})f(x),$

where $f(x)=[\sin^{-1}(1/\sqrt{x})]^{2}$ for $x>1$ and
$\tau_{X}=4m_{X}^{2}/m^{2}_{h}$ with $m_{h}$ being the SM Higgs boson mass.
The current bound is given by ATLAS:2019nkf

$\mu_{\gamma\gamma}=1.00\pm 0.06,$ (22)

where the absence of other exotic Higgs decay mode is assumed.

###  III.4  Dark matter

In this paper, we expect $y_{L^{\prime}_{\mu}}\gg
y_{L^{\prime}_{e}},y_{L^{\prime}_{\tau}}$ and $y_{S_{\mu}}\gg
y_{S_{e}},y_{S_{\tau}}$ to enhance the muon $g-2$. Under this assumption,
thermally averaged cross section to explain the relic density is expanded in
terms of relative velocity $v_{rel}^{2}\approx 0.2$ as $(\sigma v_{\rm
rel})\approx a_{\rm eff}+b_{\rm eff}v^{2}_{\rm rel}+\cdots$, where $a_{\rm
eff}$ and $b_{\rm eff}$ are respectively the coefficients of $s$-wave and
$p$-wave. In our case, the cross section is $p$-wave dominant, and the
contribution by the Yukawa terms is given by

$\displaystyle(\sigma v_{\rm rel})_{Y}$
$\displaystyle\simeq\frac{|y_{L^{\prime}_{\mu}}|^{4}}{392\pi}\left[s_{c}^{4}\frac{m_{R}^{2}}{(m_{R}^{2}+M_{1}^{2})^{2}}+c_{c}^{4}\frac{m_{R}^{2}}{(m_{R}^{2}+M_{2}^{2})^{2}}+\frac{m_{R}^{2}}{(m_{R}^{2}+M_{L^{\prime}}^{2})^{2}}\right]v^{2}_{\rm
rel}$
$\displaystyle+\frac{|y_{S_{\mu}}|^{4}}{392\pi}\left[c_{c}^{4}\frac{m_{R}^{2}}{(m_{R}^{2}+M_{1}^{2})^{2}}+s_{c}^{4}\frac{m_{R}^{2}}{(m_{R}^{2}+M_{2}^{2})^{2}}\right]v^{2}_{\rm
rel}.$ (23)

We expect the cross section from the Yukawa terms are not dominant as well as
Higgs potential (except the pole at around half of Higgs mass $\sim 63$ GeV)
due to the constraint of direct detection search. In fact, we find our cross
section reaches $10^{-10}$ GeV-2 at most, that is below the value of cross
section to explain relic density by one order magnitude through our numerical
estimation. Therefore, we expect the current relic density is given via gauge
interactions such as $2\eta_{R}\to W^{+}W^{-}/ZZ^{*}$ via
$\eta^{\pm}/\eta_{I}$. The detailed analysis has already been done by ref.
Hambye:2009pw  and they tell us the DM mass, which satisfies the constraint of
the relic abundance is about 534 GeV (in addition to the Higgs pole mass $\sim
63$ GeV).

Direct detection: We have a spin independent scattering cross section
constrained by several experiments such as XENON1T XENON:2015gkh . Since there
is a mass difference between $\eta_{R}$ and $\eta_{I}$, we can simply evade
the $Z$ portal scattering cross section so that the mass difference is taken
to be more than $100$ keV. The dominant contribution is the Higgs portal via
interactions in the scalar potential. However, we simply assume the Higgs
portal coupling $c_{\eta\eta h}$ for $\eta_{R}-\eta_{R}-h$ to be tiny enough;
$c_{\eta\eta h}\lesssim 10^{-3}$ Arcadi:2019lka , in order to evade the direct
detection via Higgs exchanges.

##  IV Numerical analysis

In this section, we perform the numerical analysis. Relevant free parameters
are randomly scanned to search for the allowed region to satisfy muon $g-2$
and oblique parameters. The ranges of the other input parameters are set to be
as follows:

$\displaystyle[y_{L^{\prime}_{\mu}},y_{S_{\mu}}]\in[0.1,1],\
[y_{L^{\prime}_{e,\tau}},y_{S_{e,\tau}}]\leq 0.001,$ $\displaystyle
m_{R}\in[70,\ 550]\ [{\rm GeV}],\ \\{M_{E,L^{\prime}}\\}\in[m_{R},\ 5000]\
[{\rm GeV}],\ $ (24) $\displaystyle m_{E}\in[10,\sqrt{4\pi}v]\ [{\rm GeV}],\
m_{I}\in[m_{R},\ \sqrt{m_{R}^{2}+\sqrt{4\pi}v^{2}}]\ [{\rm GeV}]$

where the lower bound on $m_{R}$ comes from LEP bound, we take the upper
bounds on $m_{I}$ and $m_{E}$ so that the perturbative limit for dimensionless
couplings be less than $\sqrt{4\pi}$, and we choose $M_{1}\leq M_{2}$, and
take $m_{R}\leq(M_{1},m_{I})$ so that $m_{R}$ is DM.

Figure 2: Oblique parameters $\Delta S$ and $\Delta T$ in light of the W boson
anomaly, where the blue circle represents the allowed region at $68\%$
confidence level (CL), green $90\%$ CL, and red $99\%$ CL. The black plots are
the ones satisfying the muon $g-2$ at 3$\sigma$ interval. The orange(white)
plots satisfy the lighter(heavier) mass of DM at the region of
$62.5-63(525.5-542.5)$ GeV in addition to the muon $g-2$ at 3$\sigma$
interval.

Fig. 2 shows oblique parameters $\Delta S$ and $\Delta T$ in light of the W
boson anomaly, where the blue circle represents the allowed region at $68\%$
confidence level (CL), green $90\%$ CL, and red $99\%$ CL. The black plots are
the ones satisfying the muon $g-2$ at 3$\sigma$ interval. The orange(white)
plots satisfy the lighter(heavier) mass of DM at the region of
$62.5-63(525.5-542.5)$ GeV in addition to the muon $g-2$ at 3$\sigma$
interval. It implies that the lighter DM would be in favor of the W boson
anomaly. Hereafter, we depict figures within the range of oblique parameters.

Figure 3: The left(right) figure is muon $g-2$ at $3\sigma$ interval in terms
of $m_{R}-m_{I}$($m_{R}$), where the red points show the allowed region of
oblique parameters within $68\%$ confidential level (CL), green $90\%$ CL, and
blue $99\%$ CL. The solid black line represents the BF of $\Delta a_{\mu}$,
the dotted lines $1\sigma$, and the dot-dashed lines $2\sigma$.

The left(right) of Fig. 3 shows muon $g-2$ at $3\sigma$ interval in terms of
$m_{R}-m_{I}$($m_{R}$), where the red points show the allowed region of
oblique parameters within $68\%$ CL, green $90\%$ CL, and blue $99\%$ CL. The
solid black line represents the best fit (BF) of $\Delta a_{\mu}$, the dotted
lines $1\sigma$, and the dot-dashed lines $2\sigma$. The smaller mass
difference tends to get smaller muon $g-2$, that is trivial result since the
muon $g-2$ is proportional to the mass difference. The smaller DM mass
provides a wider allowed region that is also a natural consequence of the
formula of muon $g-2$. When the DM mass and the mass difference are heavier,
the allowed space is narrower since they conflict with the oblique parameters.

Figure 4: The left(right) figure is muon $g-2$($M_{2}$) at $3\sigma$ interval
in terms of $M_{1}$, where the color and line legends are the same as in the
previous figure.

The left(right) Fig. 4 shows muon $g-2$ at $3\sigma$ interval in terms of
$m_{R}-m_{I}$($m_{R}$), where the color and line legends are the same as in
the previous figure. Similar to the previous figure, the smaller mass of
$M_{1}$ tends to have more points, and the allowed mass ranges are found to be
$100\ {\rm GeV}\lesssim M_{1}\lesssim 1700$ GeV and $100\ {\rm GeV}\lesssim
M_{1}\lesssim 5000$ GeV within our input ranges.

Figure 5: $\mu_{\gamma\gamma}$ in terms of $\Delta T$ with $\lambda_{H\eta\ll
1}$, where the colors are the same as the previous figure. The region between
the dashed line corresponds to the allowed region at 90$\%$ CL.

We show $\mu_{\gamma\gamma}$ for our allowed parameter sets as a function of
$\Delta T$ where the colors are the same as the previous figures and the
region between dashed line corresponds to allowed region at 90$\%$ CL. Here we
consider $\lambda_{H\eta}\ll 1$ to see the effect of exotic charged leptons
for illustration. The deviation from the SM ($\mu_{\gamma\gamma}=1$) can be
sizable. Note that the points outside the allowed region can be shifted inside
the allowed region by choosing the value of $\lambda_{H\eta}$ to compensate
charged fermion contribution.

Finally, we show a benchmark point to satisfy DM, muon $g-2$ at $3\sigma$
interval, and oblique parameters at $68\%$ CL arising from the W boson mass
anomaly in Tab. 2. Here we also comment on constraints from charged particle
production at the collider experiments since our $E_{1}$ is as light as $105$
GeV. At the LHC, it can be pair produced as $pp\to E_{1}^{+}E_{1}^{-}$ via
electroweak interactions. Then $E_{1}$ decays into $\eta_{R}\mu$ via Yukawa
interaction. As a constraint, we refer to the search for slepton pair
production in which slepton decays into charged lepton and neutralino(DM)
ATLAS:2019lff . We find that $105$ GeV slepton and $\sim 60$ GeV neutralino
region is still allowed, and our benchmark point is still safe from collider
constraint. Heavier charged particle has more complicated decay chains and we
need detailed analysis to obtain discovery potential at collider experiments.
A more detailed discussion of collider physics is beyond our scope of this
work and it is left as future work.

$[m_{R},m_{I},m_{\eta^{\pm}}]/{\rm GeV}$ | $[62.9,237,288]$
---|---
$[M_{E},M_{L^{\prime}},m_{E}]/{\rm GeV}$ | $[1101,124.4,137.7]$
$[M_{1},M_{2}]/{\rm GeV}$ | $[105.4,1120]$
$[y_{L^{\prime}_{\mu}},y_{S_{\mu}}]$ | $[0.0194,-0.307]$
$\Delta m^{2}_{\rm atm}$ | $2.52\times 10^{-3}{\rm eV}^{2}$
$(\sigma v_{\rm rel})_{Y}\ {\rm GeV^{2}}$ | $1.29\times 10^{-14}$
$\Delta a_{\mu}$ | $9.73\times 10^{-10}$
$\Delta T$ | $0.100$
$\Delta S$ | $0.00248$
$\mu_{\gamma\gamma}$ | $0.93$
Table 2: Numerical benchmark points to satisfy DM, muon $g-2$ at $3\sigma$
interval, and oblique parameters at $68\%$ CL.

##  V Summary and Conclusions

We have explained the W boson mass anomaly reported by CDFII collaboration
that might suggest NP. For this purpose, we have introduced exotic fermions;
one isospin doublet vector fermion, one isospin singlet singly-charged vector
fermion, and an isospin doublet inert boson. These new fermions also provide
us sizable muon $g-2$ due to no chiral suppression. Note that we need a
nonzero mass difference between $\eta_{R}$ and $\eta_{I}$; otherwise muon
$g-2$ is zero. This mass difference affects the oblique parameters, especially
$\Delta T$, whose shift from zero is the direct consequence of this W boson
mass anomaly. The doublet fermion negatively contribute to $\Delta T$ (as well
as $\Delta S$), and the doublet boson can positively contributes to it.
Because the current deviation tends to shift to the positive side of $\Delta
T$ as shown in Fig. 2, the doublet boson is more important. On the other hand,
the too big mass difference between $\eta_{R}$ and $\eta_{I}$ could conflict
with these oblique parameters. Therefore, oblique parameters and muon $g-2$
might conflict with each other, and we need to test possible phenomenological
constraints taking into account both scalar and fermion contributions.

Then, we have performed the numerical analysis to satisfy both, and we have
found that there still exists an allowed region, as shown in Figs. 3 and 4. In
addition, we have discussed a DM candidate assigned by $\eta_{R}$, and we have
found that lighter DM mass($\sim 63$ GeV) is favored through our analysis.
Here, we have considered a situation in which the thermally averaged cross
section to explain the relic density dominantly arises from gauge
interactions. In fact, we have confirmed that our cross sections from Yukawas
$y_{L^{\prime}}$ and $y_{S}$ are suppressed enough to be neglected. Under this
condition, we have obtained another result that muon $g-2$ is lower a bit that
is less than $10^{-9}$ when we take into consideration of DM.

## Acknowledgments

KIN was supported by JSPS Grant-in-Aid for Scientific Research (A) 18H03699,
(C) 21K03562, (C) 21K03583, Okayama Foundation for Science and Technology, and
Wesco Scientific Promotion Foundation. The work of H.O. is supported by an
appointment to the JRG Program at the APCTP through the Science and Technology
Promotion Fund and Lottery Fund of the Korean Government, and also by the
Korean Local Governments - Gyeongsangbuk-do Province and Pohang City H. O. is
sincerely grateful for the KIAS members. The work was supported in part by the
Fundamental Research Funds for the Central Universities (T. N.).

##  Appendix A Vector-like fermion contributions to electroweak vacuum
polarization diagrams

To estimate the oblique parameters, we evaluate the electroweak vacuum
polarization diagrams with vector-like fermions and then summarize the
analytic form of the contributions. The contributions of vacuum polarizations
for $Z$ and $W$ to the oblique parameters are given

$\displaystyle\Pi_{Z}^{\mu\nu}=g^{\mu\nu}\frac{e^{2}}{c_{W}^{2}s_{W}^{2}}(\Pi_{33}(q^{2})-2s_{W}^{2}\Pi_{3Q}(q^{2})-s_{W}^{4}\Pi_{QQ}(q^{2})),$
(25)
$\displaystyle\Pi^{\mu\nu}_{W}=g^{\mu\nu}\frac{e^{2}}{s_{W}^{2}}\Pi_{\pm}(q^{2}),$
(26)

where $q$ is four momentum carried by gauge bosons. Those with the non-zero
contribution to $\Pi_{33}(q^{2})$ are listed below:

$\displaystyle\Pi^{E_{1}E_{1}}_{33}(q^{2})=-\frac{s_{c}^{4}}{16\pi^{2}}F(q^{2},M_{1}^{2},M_{1}^{2}),$
(27)
$\displaystyle\Pi^{E_{2}E_{2}}_{33}(q^{2})=-\frac{c_{c}^{4}}{16\pi^{2}}F(q^{2},M_{2}^{2},M_{2}^{2}),$
(28)
$\displaystyle\Pi^{E_{1}E_{2}}_{33}(q^{2})=-\frac{s_{c}^{2}c_{c}^{2}}{8\pi^{2}}F(q^{2},M_{1}^{2},M_{2}^{2}),$
(29)
$\displaystyle\Pi^{N^{\prime}N^{\prime}}_{33}(q^{2})=-\frac{1}{16\pi^{2}}F(q^{2},M_{N^{\prime}}^{2},M_{N^{\prime}}^{2}),$
(30)

where the superscripts on the left side represent the particles inside the
vacuum polarization diagrams. $F(q^{2},m^{2},m^{\prime 2})$ is the loop
function given by

$\displaystyle F(q^{2},m^{2},m^{\prime 2})=$
$\displaystyle\int_{0}^{1}dxdy\delta(1-x-y)\left(\frac{1}{\epsilon_{\rm\overline{MS}}}-\ln\left(\frac{\Delta}{\mu^{2}}\right)\right)$
$\displaystyle\times(2x(1-x)q^{2}-xm^{2}-ym^{\prime 2}+mm^{\prime}),$ (31)
$\displaystyle\Delta=$ $\displaystyle-x(1-x)q^{2}+xm^{2}+ym^{\prime
2},\quad\frac{1}{\epsilon_{\rm\overline{MS}}}\equiv\frac{2}{\epsilon}-\gamma-\ln(4\pi),$
(32)

where $\mu$ is auxiliary parameter with mass dimension. Once the
$S,T$-parameters are computed from Eqs. (12) and (13) in Sec. III.1, the $\mu$
dependence cancels out and thus does not remain.

In a similar way, one finds a non-zero contribution for $\Pi_{3Q,QQ,\pm}$:

$\displaystyle\Pi^{E_{1}E_{1}}_{3Q}(q^{2})=-\frac{s_{c}^{2}}{8\pi^{2}}F(q^{2},M_{1}^{2},M_{1}^{2}),$
(33)
$\displaystyle\Pi^{E_{2}E_{2}}_{3Q}(q^{2})=-\frac{c_{c}^{2}}{8\pi^{2}}F(q^{2},M_{2}^{2},M_{2}^{2}),$
(34)
$\displaystyle\Pi^{E_{1}E_{1}}_{QQ}(q^{2})=-\frac{1}{4\pi^{2}}F(q^{2},M_{1}^{2},M_{1}^{2}),$
(35)
$\displaystyle\Pi^{E_{2}E_{2}}_{QQ}(q^{2})=-\frac{1}{4\pi^{2}}F(q^{2},M_{2}^{2},M_{2}^{2}),$
(36)
$\displaystyle\Pi^{E_{1}N^{\prime}}_{\pm}(q^{2})=-\frac{s_{c}^{2}}{8\pi^{2}}F(q^{2},M_{1}^{2},M_{N^{\prime}}^{2}),$
(37)
$\displaystyle\Pi^{E_{2}N^{\prime}}_{\pm}(q^{2})=-\frac{c_{c}^{2}}{8\pi^{2}}F(q^{2},M_{2}^{2},M_{N^{\prime}}^{2}).$
(38)

Then the $S$\- and $T$-parameters can be obtained from Eqs. (12) and (13) in
Sec. III.1 and those equations show that the divergent part proportional to
$\epsilon_{\rm SM}$ is canceled if the oblique parameters are computed.

## References

  * (1) T. Aaltonen et al. [CDF], Science 376 (2022) no.6589, 170-176 doi:10.1126/science.abk1781
  * (2) M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98 (2018) no.3, 030001 doi:10.1103/PhysRevD.98.030001
  * (3) Y. Z. Fan, T. P. Tang, Y. L. S. Tsai and L. Wu, [arXiv:2204.03693 [hep-ph]].
  * (4) C. T. Lu, L. Wu, Y. Wu and B. Zhu, [arXiv:2204.03796 [hep-ph]].
  * (5) P. Athron, A. Fowlie, C. T. Lu, L. Wu, Y. Wu and B. Zhu, [arXiv:2204.03996 [hep-ph]].
  * (6) G. W. Yuan, L. Zu, L. Feng and Y. F. Cai, [arXiv:2204.04183 [hep-ph]].
  * (7) A. Strumia, [arXiv:2204.04191 [hep-ph]].
  * (8) J. M. Yang and Y. Zhang, [arXiv:2204.04202 [hep-ph]].
  * (9) J. de Blas, M. Pierini, L. Reina and L. Silvestrini, [arXiv:2204.04204 [hep-ph]].
  * (10) G. G. di Cortona and E. Nardi, [arXiv:2204.04227 [hep-ph]].
  * (11) X. K. Du, Z. Li, F. Wang and Y. K. Zhang, [arXiv:2204.04286 [hep-ph]].
  * (12) T. P. Tang, M. Abdughani, L. Feng, Y. L. S. Tsai and Y. Z. Fan, [arXiv:2204.04356 [hep-ph]].
  * (13) G. Cacciapaglia and F. Sannino, [arXiv:2204.04514 [hep-ph]].
  * (14) M. Blennow, P. Coloma, E. Fernández-Martínez and M. González-López, [arXiv:2204.04559 [hep-ph]].
  * (15) K. Sakurai, F. Takahashi and W. Yin, [arXiv:2204.04770 [hep-ph]].
  * (16) J. Fan, L. Li, T. Liu and K. F. Lyu, [arXiv:2204.04805 [hep-ph]].
  * (17) X. Liu, S. Y. Guo, B. Zhu and Y. Li, [arXiv:2204.04834 [hep-ph]].
  * (18) H. M. Lee and K. Yamashita, [arXiv:2204.05024 [hep-ph]].
  * (19) Y. Cheng, X. G. He, Z. L. Huang and M. W. Li, [arXiv:2204.05031 [hep-ph]].
  * (20) H. Song, W. Su and M. Zhang, [arXiv:2204.05085 [hep-ph]].
  * (21) E. Bagnaschi, J. Ellis, M. Madigan, K. Mimasu, V. Sanz and T. You, [arXiv:2204.05260 [hep-ph]].
  * (22) A. Paul and M. Valli, [arXiv:2204.05267 [hep-ph]].
  * (23) H. Bahl, J. Braathen and G. Weiglein, [arXiv:2204.05269 [hep-ph]].
  * (24) P. Asadi, C. Cesarotti, K. Fraser, S. Homiller and A. Parikh, [arXiv:2204.05283 [hep-ph]].
  * (25) L. Di Luzio, R. Gröber and P. Paradisi, [arXiv:2204.05284 [hep-ph]].
  * (26) P. Athron, M. Bach, D. H. J. Jacob, W. Kotlarski, D. Stöckinger and A. Voigt, [arXiv:2204.05285 [hep-ph]].
  * (27) J. Gu, Z. Liu, T. Ma and J. Shu, [arXiv:2204.05296 [hep-ph]].
  * (28) J. J. Heckman, [arXiv:2204.05302 [hep-ph]].
  * (29) K. S. Babu, S. Jana and V. P. K., [arXiv:2204.05303 [hep-ph]].
  * (30) Y. Heo, D. W. Jung and J. S. Lee, [arXiv:2204.05728 [hep-ph]].
  * (31) X. K. Du, Z. Li, F. Wang and Y. K. Zhang, [arXiv:2204.05760 [hep-ph]].
  * (32) K. Cheung, W. Y. Keung and P. Y. Tseng, [arXiv:2204.05942 [hep-ph]].
  * (33) A. Crivellin, M. Kirk, T. Kitahara and F. Mescia, [arXiv:2204.05962 [hep-ph]].
  * (34) M. Endo and S. Mishima, [arXiv:2204.05965 [hep-ph]].
  * (35) T. Biekötter, S. Heinemeyer and G. Weiglein, [arXiv:2204.05975 [hep-ph]].
  * (36) R. Balkin, E. Madge, T. Menzo, G. Perez, Y. Soreq and J. Zupan, [arXiv:2204.05992 [hep-ph]].
  * (37) N. V. Krasnikov, [arXiv:2204.06327 [hep-ph]].
  * (38) Y. H. Ahn, S. K. Kang and R. Ramos, [arXiv:2204.06485 [hep-ph]].
  * (39) X. F. Han, F. Wang, L. Wang, J. M. Yang and Y. Zhang, [arXiv:2204.06505 [hep-ph]].
  * (40) M. D. Zheng, F. Z. Chen and H. H. Zhang, [arXiv:2204.06541 [hep-ph]].
  * (41) J. Kawamura, S. Okawa and Y. Omura, [arXiv:2204.07022 [hep-ph]].
  * (42) A. Ghoshal, N. Okada, S. Okada, D. Raut, Q. Shafi and A. Thapa, [arXiv:2204.07138 [hep-ph]].
  * (43) P. F. Perez, H. H. Patel and A. D. Plascencia, [arXiv:2204.07144 [hep-ph]].
  * (44) M. E. Peskin and T. Takeuchi, Phys. Rev. D 46 (1992), 381-409 doi:10.1103/PhysRevD.46.381
  * (45) M. E. Peskin and T. Takeuchi, Phys. Rev. Lett. 65 (1990), 964-967 doi:10.1103/PhysRevLett.65.964
  * (46) P. Ko, T. Nomura and H. Okada, [arXiv:2110.10513 [hep-ph]].
  * (47) C. W. Chiang and H. Okada, Int. J. Mod. Phys. A 34 (2019) no.20, 1950106 doi:10.1142/S0217751X19501069 [arXiv:1711.07365 [hep-ph]].
  * (48) R. Barbieri, L. J. Hall and V. S. Rychkov, Phys. Rev. D 74 (2006), 015007 doi:10.1103/PhysRevD.74.015007 [arXiv:hep-ph/0603188 [hep-ph]].
  * (49) B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126 (2021) no.14, 141801 doi:10.1103/PhysRevLett.126.141801 [arXiv:2104.03281 [hep-ex]].
  * (50) T. Hambye, F. S. Ling, L. Lopez Honorez and J. Rocher, JHEP 07 (2009), 090 [erratum: JHEP 05 (2010), 066] doi:10.1007/JHEP05(2010)066 [arXiv:0903.4010 [hep-ph]].
  * (51) E. Aprile et al. [XENON], JCAP 04 (2016), 027 doi:10.1088/1475-7516/2016/04/027 [arXiv:1512.07501 [physics.ins-det]].
  * (52) G. Arcadi, A. Djouadi and M. Raidal, Phys. Rept. 842 (2020), 1-180 doi:10.1016/j.physrep.2019.11.003 [arXiv:1903.03616 [hep-ph]].
  * (53) G. Aad et al. [ATLAS], Phys. Rev. D 101 (2020) no.1, 012002 doi:10.1103/PhysRevD.101.012002 [arXiv:1909.02845 [hep-ex]].
  * (54) G. Aad et al. [ATLAS], Eur. Phys. J. C 80 (2020) no.2, 123 doi:10.1140/epjc/s10052-019-7594-6 [arXiv:1908.08215 [hep-ex]].


# Spatial Censored Regression Models in R: The CensSpatial Package

José Alejandro Ordoñez
Department of Statistics
State University of Campinas
Campinas, Brazil 13083-970
<EMAIL_ADDRESS>
&Christian E. Galarza
Department of Mathematics
Escuela Superior Politécnica del Litoral
Guayaquil, Ecuador 090112
<EMAIL_ADDRESS>
&Victor H. Lachos
Department of Statistics
University of Connecticut
Connecticut, USA 06269
<EMAIL_ADDRESS>

###### Abstract

CensSpatial is an R package for analyzing spatial censored data through linear
models. It offers a set of tools for simulating, estimating, making
predictions, and performing local influence diagnostics for outlier detection.
The package provides four algorithms for estimation and prediction. One of
them is based on the stochastic approximation of the EM (SAEM) algorithm,
which allows easy and fast estimation of the parameters of linear spatial
models when censoring is present. The package provides worthy measures to
perform diagnostic analysis using the Hessian matrix of the completed log-
likelihood function. This work is divided into two parts. The first part
discusses and illustrates the utilities that the package offers for estimating
and predicting spatial censored data. The second one describes the valuable
tools to perform diagnostic analysis. Several examples in spatial
environmental data are also provided.

##  1 Introduction

Spatial data can be described as observations that are spatially located in a
coordinated area. Their analysis constitutes an area of statistics that is
characterized by the spatial correlation among individuals where the interest
is to predict the response at non-sampling sites. Usually, spatial analysis
assumes that observations are fully observed, but this is not always possible.
Sometimes data are subject to upper or lower (or both) detection limits,
beyond which they are not quantifiable (Schelin and Sjöstedt-de Luna, 2014).
An example of this can be seen in environmental (spatial) monitoring of
different variables which often involves left-censored observations falling
below the minimum limit of detection (LOD) of the instruments used to quantify
them.

Methods for handling this type of data have been widely discussed in the
statistical literature. Rathbun (2006) applied a Robbins-Monro stochastic
approximation to estimate the parameters of a spatial regression model, obtain
conditional simulations of left-censored observations and then obtain a
predictor for data at unsampled sites by taking the weighted mean of kriging
predictors via importance sampling. Schelin and Sjöstedt-de Luna (2014)
proposed a semi-naive method that determines imputed values at censored
locations in an iterative algorithm together with variogram estimation. They
compared its predictive performance through a simulation study with Rathbun’s,
Naive 1 and Naive 2 methods, which involve replacing the censored data by the
LOD or the LOD/2 respectively. Militino and Ugarte (1999) developed an EM-type
algorithm for maximum likelihood (ML) estimation in censored spatial data.
However, this approach suffers from several drawbacks that restrict its
applicability, such as the computation of non-closed form expressions, which
depend on high dimensional integrals.

The EM algorithm proposed by Dempster et al. (1977) is an iterative procedure
for finding the maximum likelihood (ML) estimators of models that involves
unobserved latent variables. In our context, an EM-type algorithm uses the
observed and censored data to obtain estimates from the mean and correlation
structures of a spatial model. In some cases, this algorithm cannot be carried
out because the calculation of the expectations involved in E-step cannot be
computed analytically. An alternative is to approximate those expectations by
Monte Carlo methods to replace the traditional E-step. This procedure is the
so-called Monte Carlo Markov chain (MCEM) algorithm. Even though it is
possible to approximate these expectations, this can involve heavy
computational effort when some generating methods are used to sample from the
conditional latent variable distributions at each iteration.

An alternative to deal with this problem is the SAEM algorithm proposed by
Delyon et al. (1999). This algorithm replaces the E-step by a stochastic
approximation obtained using simulated data, while the M-step remains
unchanged. In the framework of spatial models, Jank (2006) showed that the
computational effort of SAEM is much smaller and it reaches convergence in
just a fraction of the simulation iterations when compared with the MCEM
algorithm. This is due to the fact that a memory effect is contained in the
SAEM method, in which the previous simulations are considered in the
computation of the posterior ones.

Another important topic in spatial censored data analysis is the study of
influential observations. Influence analysis is an important statistical tool
because it can assess the goodness-of-fit and influential observations that
can distort the parameter estimates leading in some cases to erroneous
inference (Giménez and Galea, 2013). For influential observation detection, we
discuss two approaches; the first one is the case-deletion approach (Cook,
1977), an intuitively appealing method which has been applied to many
statistical models. The second one is local influence diagnostics (Cook,
1986). This technique assesses the stability of the estimation outputs with
respect to the model inputs and has recently received attention in the
literature on spatial models (see, e.g., De Bastiani et al., 2014). Zhu and
Lee (2001) developed an approach to perform local influence analysis for
general statistical models with missing data, based on the 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 the Q-displacement function, following the
approach of Zhu et al. (2001). So, we develop here methods to obtain case-
deletion measures and local influence measures by using the method of Zhu et
al. (2001) (see also Lee and Xu, 2004; Zhu and Lee, 2001) in the context of
spatial censored models.

The aim of this paper is to present the utilities of the R package
CensSpatial. The rest of the paper is organized as follows. In Section 2, we
derive all the theoretical framework for the estimation and prediction
functions for spatial censored data implemented in the proposed package, along
with the R functions for calling these methods. Also we present some useful
graphical tools for interpretation as well as an application to real data.
Section 4 presents the local influence diagnostic methods to detect
influential observations. Graphical tools and an application are also provided
in this section. Further details such as first and second derivatives of the
spatial correlation, can be found in the Appendix.

##  2 Estimation and prediction utilities

First, we specify the spatial censored linear (SCL) model in a general sense
and then we describe all the estimation and prediction tools implemented in
the CensSpatial package.

###  2.1 Background of spatial linear models for censored response

As in De Bastiani et al. (2014), we will consider a Gaussian model with a
linear specification for the spatial trend, which allows the inclusion of
polynomial trends or more generally, spatially referenced covariates. The
linear model is defined as

$\mathbf{Z}=\mathbf{X}\mbox{\boldmath$\beta$}+\mbox{\boldmath$\epsilon$},$ (1)

where $\mathbf{X}$ and $\beta$ represent the trend matrix and the parameters
considered for the mean structure respectively, and
$\mbox{\boldmath$\epsilon$}\sim N_{n}(\mathbf{0},\mbox{\boldmath$\Sigma$})$ is
the stochastic component. We consider for the spatial case, $\mathbf{X}$ is a
full rank matrix and $\Sigma$ is a non-singular matrix with the form
$\mbox{\boldmath$\Sigma$}=\sigma^{2}\mbox{\boldmath$\Psi$}$, with
$\mbox{\boldmath$\Psi$}=(\nu^{2}\mathbf{I}_{n}+\mathbf{R}(\phi))$, where
$\mathbf{I}_{n}$ represents an $n\times n$ identity matrix, $\mathbf{R}(\phi)$
is the correlation matrix for the error term and $\nu^{2}=\tau^{2}/\sigma^{2}$
(see Appendix). We also assume that the response $Z_{i}$ is not fully observed
for all area $i$, so we call $(V_{i},C_{i})$ the observed data for the $i$th
area, while $V_{i}$ represents either an uncensored observation
$(V_{i}=V_{0i})$ or the LOD of censoring level $(V_{i}=(V_{1i},V_{2i}))$ and
$C_{i}$ is the censoring indicator such that:

$\displaystyle C_{i}=\left\\{\begin{array}[]{ccc}1&if&V_{1i}\leq Z_{i}\leq
V_{2i}\,,\\\ 0&if&Z_{i}=V_{0i}\,.\end{array}\right.$ (4)

Note that when $C_{i}=1$, if $Z_{i}\in(-\infty,V_{2i}]$, then we get a left
censored SCL model (Toscas, 2010), and if $Z_{i}\in[V_{1i},\infty)$, then we
get a right censored SCL model. The model defined in (1)-(4) is called the
spatial censored linear (SCL) model.

##  3 Likelihood and SAEM implementation

To compute the likelihood function associated with the SCL model, the observed
and censored components of $\mathbf{Z}$ must be treated separately. Let
$\mathbf{Z}^{o}$ be the $n^{o}$-vector of observed responses and
$\mathbf{Z}^{c}$ be the $n^{c}$-vector of censored observations, with
$(n=n^{o}+n^{c})$, such that $C_{i}=0$ for all elements in $\mathbf{Z}^{o}$,
and $C_{i}=1$ for all elements in $\mathbf{Z}^{c}$. After reordering,
$\mathbf{Z}$, V, $\mathbf{X}$, and $\Sigma$ can be partitioned as follows:

$\mathbf{Z}=vec(\mathbf{Z}^{o},\mathbf{Z}^{c}),\quad\textbf{V}=vec(\mathbf{V}^{o},\mathbf{V}^{c}),\quad\mathbf{X}^{\top}=[\mathbf{X}^{o\top},\mathbf{X}^{c\top}]\quad{\rm
and}\quad\mbox{\boldmath$\Sigma$}=\left[\begin{array}[]{cc}\mbox{\boldmath$\Sigma$}^{oo}&\mbox{\boldmath$\Sigma$}^{oc}\\\
\mbox{\boldmath$\Sigma$}^{co}&\mbox{\boldmath$\Sigma$}^{cc}\end{array}\right],$

where $vec(\cdot)$ denotes the function which stacks vectors or matrices
having the same number of columns. Consequently, $\mathbf{Z}^{o}\sim
N_{n^{o}}(\mathbf{X}^{o}\mbox{\boldmath$\beta$},\mbox{\boldmath$\Sigma$}^{oo})$,
$\mathbf{Z}^{c}|\mathbf{Z}^{o}\sim
N_{n^{c}}(\mbox{\boldmath$\mu$},\mathbf{S}),$ where
$\mbox{\boldmath$\mu$}=\mathbf{X}^{c}\mbox{\boldmath$\beta$}+\mbox{\boldmath$\Sigma$}^{co}(\mbox{\boldmath$\Sigma$}^{oo})^{-1}(\mathbf{Z}^{o}-\mathbf{X}^{o}\mbox{\boldmath$\beta$})$
and
$\mathbf{S}=\mbox{\boldmath$\Sigma$}^{cc}-\mbox{\boldmath$\Sigma$}^{co}(\mbox{\boldmath$\Sigma$}^{oo})^{-1}\mbox{\boldmath$\Sigma$}^{oc}$.
Now, let $\phi_{n}(\mathbf{u};\mathbf{a},\mathbf{A})$ be the _pdf_ of
$N_{n}(\mathbf{a},\mathbf{A})$ evaluated at $\mathbf{u}$. From Vaida and Liu
(2009) and Jacqmin-Gadda et al. (2000), the likelihood function (using
conditional probability arguments) is given by

$\displaystyle
L(\mbox{\boldmath$\theta$})=\phi_{n^{o}}(\mathbf{Z}^{o};\mathbf{X}^{o}\mbox{\boldmath$\beta$},\mbox{\boldmath$\Sigma$}^{oo})P(\mathbf{Z}^{c}\in\mathbb{V}^{c}|\mathbf{Z}^{o}),$
(5)

where

$\mathbb{V}^{c}=\\{\mathbf{Z}^{c}=(Z_{1}^{c},\ldots,Z_{n^{c}}^{c})^{\top}|V_{11}^{c}\leq
Z_{1}^{c}\leq V_{21}^{c},\ldots,V_{1n^{c}}^{c}\leq Z_{n^{c}}^{c}\leq
V_{2n^{c}}^{c}\\}$

and $P(\mathbf{u}\in\mathbb{A}|\mathbf{Z}^{o})$ denotes the conditional
probability of $\mathbf{u}$ being in the set $\mathbb{A}$ given the observed
response. This function can be evaluated without much computational burden
through the mvtnorm routine available in the namesake R package (see Genz et
al., 2020). Now, we can compare models using a likelihood based criterion,
such as the Akaike (AIC; Akaike, 1974) and Schwarz (BIC; Schwarz, 1978)
information criteria.

###  3.1 SAEM algorithm for censored spatial data

We propose an EM-type (SAEM) algorithm by considering $\mathbf{Z}$ as a
missing data or latent variable. In the estimation, the reading $\mathbf{V}$
at the censored observations ($C_{i}=1$) is treated as hypothetical missing
data, and augmented with the observed dataset
$\mathbf{Z}_{c}=(\mathbf{C}^{\top},\mathbf{V}^{\top},\mathbf{Z}^{\top})^{\top}$.
Hence, the complete-data log-likelihood function is given by:

$\displaystyle\ell_{c}(\mbox{\boldmath$\theta$})$ $\displaystyle\propto$
$\displaystyle-\frac{1}{2}\left[\log(|\mbox{\boldmath$\Sigma$}|)+(\mathbf{Z}-\mathbf{X}\mbox{\boldmath$\beta$})^{\top}\mbox{\boldmath$\Sigma$}^{-1}(\mathbf{Z}-\mathbf{X}\mbox{\boldmath$\beta$})\right]+cte,$
(6)

with $cte$ being a constant, independent of the parameter vector $\theta$.
Given the current estimate
$\mbox{\boldmath$\theta$}=\widehat{\mbox{\boldmath$\theta$}}^{(k)}$, the
E-step computes the conditional expectation of the complete data log-
likelihood function, i.e.,
$Q(\mbox{\boldmath$\theta$}|\widehat{\mbox{\boldmath$\theta$}}^{(k)})=E[\ell_{c}(\mbox{\boldmath$\theta$}|\mathbf{Z}_{c})|\mathbf{V},\mathbf{C},\widehat{\mbox{\boldmath$\theta$}}^{(k)}]$.
Denote the two conditional first moments for the response $\mathbf{Z}$ as
$\widehat{\mathbf{Z}}^{(k)}=E\\{\displaystyle\mathbf{Z}|\textbf{V},\mathbf{C},\widehat{\mbox{\boldmath$\theta$}}^{(k)}\\}$
and
$\widehat{\mathbf{Z}\mathbf{Z}^{\top}}^{(k)}=E\\{\displaystyle\mathbf{Z}\mathbf{Z}^{\top}|\textbf{V},\mathbf{C},\widehat{\mbox{\boldmath$\theta$}}^{(k)}\\}$.
For the SCL model, we have

$\displaystyle
Q(\mbox{\boldmath$\theta$}|\widehat{\mbox{\boldmath$\theta$}}^{(k)})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\left[\log(|\mbox{\boldmath$\Sigma$}|)+\widehat{A}^{(k)}\right],$

where

$\widehat{A}^{(k)}=tr\big{(}\widehat{\mathbf{Z}\mathbf{Z}^{\top}}^{(k)}\mbox{\boldmath$\Sigma$}^{-1}\big{)}-2\widehat{\mathbf{Z}}^{(k)\top}\mbox{\boldmath$\Sigma$}^{-1}\mathbf{X}\mbox{\boldmath$\beta$}+\mbox{\boldmath$\beta$}^{\top}\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}^{-1}\mathbf{X}\mbox{\boldmath$\beta$}.$

$\mathbf{Z}$ is not fully observed, so the components of
$\widehat{\mathbf{Z}}^{(k)}$ and $\widehat{\mathbf{Z}\mathbf{Z}}^{(k)}$,
corresponding to $C_{i}=1$ will be estimated by the first two moments of the
truncated normal distribution respectively. When $C_{i}=0$, these components
can be obtained directly from the observed values, i.e,
$\widehat{\mathbf{Z}}^{(k)}=\mathbf{Z}^{o}$ and
$\widehat{\mathbf{Z}\mathbf{Z}}^{(k)}=\mathbf{Z}^{o}\mathbf{Z}^{o\top}$.
Although these expectations exhibit closed forms (as functions of multinormal
probabilities; for further details, please refer to Arismendi (2013)), the
calculation is computationally expensive requiring high-dimensional numerical
integrations, resulting in convergence issues when the proportion of censored
observations is non-negligible. At each iteration, the SAEM algorithm
successively simulates from the conditional distribution of the latent
variable, and updates the unknown parameters of the model. Thus, at iteration
$k$ , the SAEM proceeds as follows:

  * •

Step E-1 (sampling): Sample $\mathbf{Z}^{c}$ from a truncated normal
distribution, denoted by
$TN_{n^{c}}(\mbox{\boldmath$\mu$},\mathbf{S};\mathbb{A}^{c})$, with
$\mathbb{A}^{c}=\\{\mathbf{Z}^{c}=(Z_{1}^{c},\ldots,Z_{n^{c}}^{c})^{\top}|V_{11}^{c}\leq
Z_{1}^{c}\leq V_{21}^{c},\ldots,V_{1n^{c}}^{c}\leq Z_{n^{c}}^{c}\leq
V_{2n^{c}}^{c}\\},$
$\mbox{\boldmath$\mu$}=\mathbf{X}^{c}\mbox{\boldmath$\beta$}+\mbox{\boldmath$\Sigma$}^{co}(\mbox{\boldmath$\Sigma$}^{oo})^{-1}(\mathbf{Z}^{o}-\mathbf{X}^{o}\mbox{\boldmath$\beta$})$
and
$\mathbf{S}=\mbox{\boldmath$\Sigma$}^{cc}-\mbox{\boldmath$\Sigma$}^{co}(\mbox{\boldmath$\Sigma$}^{oo})^{-1}\mbox{\boldmath$\Sigma$}^{oc}$.
Here $TN_{n}(.;\mathbb{A})$ denotes the $n$-variate truncated normal
distribution in the interval $\mathbb{A}$, where
$\mathbb{A}=A_{1}\times\ldots\times A_{n}$. The new observation
$\mathbf{Z}^{(k,l)}=(Z^{(k,l)}_{1},\ldots,Z^{(k,l)}_{n^{c}},Z_{n^{c}_{i}+1},\ldots,Z_{n})$
is a sample generated for the $n^{c}$ censored cases and the observed values
(uncensored cases), for $l=1,\ldots,M.$

  * •

Step E-2 (stochastic approximation): Since we have the sequence
$\mathbf{Z}^{(k,l)}$, at the $k$-th iteration we replace the conditional
expectations $\widehat{\mathbf{Z}}^{(k)}$ and
$\widehat{\mathbf{Z}\mathbf{Z}^{\top}}^{(k)}$ by the following stochastic
approximations:

$\displaystyle\widehat{\mathbf{Z}}^{(k)}$ $\displaystyle=$
$\displaystyle\widehat{\mathbf{Z}}^{(k-1)}+\delta_{k}\left[\frac{1}{M}\sum_{l=1}^{M}\mathbf{Z}^{(k,l)}-\widehat{\mathbf{Z}}^{(k-1)}\right],$
(7) $\displaystyle\widehat{\mathbf{Z}\mathbf{Z}^{\top}}^{(k)}$
$\displaystyle=$
$\displaystyle\widehat{\mathbf{Z}\mathbf{Z}^{\top}}^{(k-1)}+\delta_{k}\left[\frac{1}{M}\sum_{l=1}^{M}\mathbf{Z}^{(k,l)}\mathbf{Z}^{(k,l)\top}-\widehat{\mathbf{Z}\mathbf{Z}^{\top}}^{(k-1)}\right],$
(8)

where $\delta_{k}$ is a smoothness parameter, i.e., a decreasing sequence of
positive numbers as defined in Kuhn and Lavielle (2004), such that
$\sum_{k=1}^{\infty}\delta_{k}=\infty$ and
$\sum_{k=1}^{\infty}\delta_{k}^{2}<\infty$. For the SAEM, the E-Step coincides
with the MCEM algorithm at the expense of a significantly smaller number of
simulations $M$ (suggested to be $M\leq 20$). For more details, see Delyon et
al. (1999).

Finally, the conditional maximization (CM) step maximizes
$Q(\mbox{\boldmath$\theta$}|\widehat{\mbox{\boldmath$\theta$}}^{(k)})$ with
respect to $\theta$ and obtains a new estimate
$\widehat{\mbox{\boldmath$\theta$}}^{(k+1)}$, as follows:

CM-Step (conditional maximization):

$\displaystyle\widehat{\mbox{\boldmath$\beta$}}^{(k+1)}$ $\displaystyle=$
$\displaystyle\big{(}\mathbf{X}^{\top}\widehat{\mbox{\boldmath$\Sigma$}}^{-1(k)}\mathbf{X}\big{)}^{-1}\mathbf{X}^{\top}\widehat{\mbox{\boldmath$\Sigma$}}^{-1(k)}\widehat{\mathbf{Z}}^{(k)},$
$\displaystyle\widehat{\sigma^{2}}^{(k+1)}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\Big{[}tr\big{(}\widehat{\mathbf{Z}\mathbf{Z}^{\top}}^{(k)}\widehat{\mbox{\boldmath$\Psi$}}^{-1(k)}\big{)}-2\widehat{\mathbf{Z}^{\top}}^{(k)}\widehat{\Psi}^{-1(k)}\mathbf{X}\widehat{\mbox{\boldmath$\beta$}}^{(k+1)}+\widehat{\mbox{\boldmath$\beta$}}^{\top(k+1)}\mathbf{X}^{\top}\widehat{\mbox{\boldmath$\Psi$}}^{-1(k)}\mathbf{X}\widehat{\mbox{\boldmath$\beta$}}^{(k+1)}\Big{]}\,\,$
$\displaystyle\widehat{\mbox{\boldmath$\alpha$}}^{(k+1)}$ $\displaystyle=$
$\displaystyle\underset{\scriptscriptstyle\mbox{\boldmath$\alpha$}\in\mathbb{R}^{+}\times\mathbb{R}^{+}}{\mathrm{argmax}}\left(-\frac{1}{2}\log(|\mbox{\boldmath$\Sigma$}|)-\frac{1}{2}\left[tr\big{(}\widehat{\mathbf{Z}\mathbf{Z}^{\top}}^{(k)}{\mbox{\boldmath$\Sigma$}}^{-1}\big{)}-2\widehat{\mathbf{Z}^{\top}}^{(k)}{\mbox{\boldmath$\Sigma$}}^{-1}\mathbf{X}\widehat{\mbox{\boldmath$\beta$}}^{(k+1)}\right.\right.$
(9) $\displaystyle+$
$\displaystyle\left.\left.\widehat{\mbox{\boldmath$\beta$}}^{\top(k+1)}\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}^{-1}\mathbf{X}\widehat{\mbox{\boldmath$\beta$}}^{(k+1)}\right]\right),$

with $\mbox{\boldmath$\alpha$}=(\nu^{2},\phi)^{\top}$. Note that
$\widehat{\tau}^{2}$ can be recovered via
$\widehat{\tau}^{2(k+1)}=\widehat{\nu}^{2(k+1)}\widehat{\sigma}^{2(k+1)}$. The
CM-step (9) can be easily accomplished using the optim routine in R. This
process is iterated until some absolute distance between two successive
evaluations of the actual log-likelihood $\ell(\mbox{\boldmath$\theta$})$,
such as
$|\ell(\widehat{\mbox{\boldmath$\theta$}}^{(k+1)})-\ell(\widehat{\mbox{\boldmath$\theta$}}^{(k)})|$
or
$|\ell(\widehat{\mbox{\boldmath$\theta$}}^{(k+1)})/\ell(\widehat{\mbox{\boldmath$\theta$}}^{(k)})-1|$,
becomes small enough.

###  3.2 Prediction

Following Diggle and Ribeiro (2007), let $\mathbf{Z}_{obs}$ denote a vector of
random variables with observed realized values, and let $\mathbf{Z}_{pred}$
denote another random variable whose realized values we would like to predict
from the observed values of $\mathbf{Z}_{obs}$. A prediction for
$\mathbf{Z}_{pred}$ could be any function of $\mathbf{Z}_{obs}$, which we
denote by $\hat{\mathbf{Z}}_{pred}=t(\mathbf{Z}_{obs})$. The mean square error
(MSE) of $\hat{\mathbf{Z}}_{pred}$ is given by
$MSE(\hat{\mathbf{Z}}_{pred})=E((\hat{\mathbf{Z}}_{pred}-\mathbf{Z}_{pred})^{2})$;
so we can use $\hat{\mathbf{Z}}_{pred}=E(\mathbf{Z}_{pred}|\mathbf{Z}_{obs})$
to find the predictor that minimizes $MSE(\mathbf{Z}_{pred})$.

For a Gaussian process, we write
$\mathbf{Z}_{pred}=(Z_{pred}(\mathbf{s}_{1}),\ldots,Z_{pred}(\mathbf{s}_{n}))$
for the unobserved values of the signal at the sampling locations
$\mathbf{s}_{1},\ldots,\mathbf{s}_{n}$. We want to predict the signal value at
an arbitrary location; thus our target for prediction is
$Z_{pred}(\mathbf{s})$. Given that $(\mathbf{Z}_{obs},\mathbf{Z}_{pred})$ is
also a multivariate Gaussian process, we can use the results in Subsection 3
to minimize the $MSE(\hat{\mathbf{Z}}_{pred})$, i.e., if
$\mathbf{X}^{*}=(\mathbf{X}_{obs},\mathbf{X}_{pred})$ is the
$(n_{obs}+n_{pred})\times p$ design matrix corresponding to
$\mathbf{Z}^{*}=(\mathbf{Z}_{obs}^{\top},\mathbf{Z}_{pred}^{\top})$, then
$\mathbf{Z}^{*}\sim
N_{n_{obs}+n_{pred}}\left(\mathbf{X}^{*}\mbox{\boldmath$\beta$},\mbox{\boldmath$\Sigma$}\right)$,
where

$\mbox{\boldmath$\Sigma$}=\begin{bmatrix}\mbox{\boldmath$\Sigma$}_{obs,obs}&\mbox{\boldmath$\Sigma$}_{obs,pred}\\\
\mbox{\boldmath$\Sigma$}_{pred,obs}&\mbox{\boldmath$\Sigma$}_{pred,pred}\end{bmatrix}$

and $\mathbf{Z}_{pred}|\mathbf{Z}_{obs}\sim
N_{n_{pred}}(\mbox{\boldmath$\mu$}_{p},\mbox{\boldmath$\Sigma$}_{p})$, where

$\mbox{\boldmath$\mu$}_{p}=\mathbf{X}_{pred}\mbox{\boldmath$\beta$}+\mbox{\boldmath$\Sigma$}_{obs,pred}(\mbox{\boldmath$\Sigma$}_{obs,obs})^{-1}(\mathbf{Z}_{obs}-\mathbf{X}_{obs}\mbox{\boldmath$\beta$})$
(10)

and

$\mbox{\boldmath$\Sigma$}_{p}=\mbox{\boldmath$\Sigma$}_{pred,pred}-\mbox{\boldmath$\Sigma$}_{pred,obs}(\mbox{\boldmath$\Sigma$}_{obs,obs})^{-1}\mbox{\boldmath$\Sigma$}_{obs,pred}.$

Then, the predictor that minimizes the $MSE$ of prediction will be the
conditional expectation in (10).

####  3.2.1 Prediction methods

Using the previous results, we will now describe the algorithm for each of the
three prediction methods that are implemented in this paper.

  1. (a)

Naive 1 and Naive 2 algorithms

We proceed as follows:

    1. 1.

Impute the censored observations using $LOD$ (called Naive 1 algorithm), or
$LOD/2$ (called Naive 2 algorithm), depending on the type of censoring (left
or right).

    2. 2.

Compute the least squares or likelihood estimates for the mean (linear trend
in our case) and covariance structure using the imputed data.

    3. 3.

Evaluate the estimate $\hat{\mbox{\boldmath$\theta$}}$ from step 2 in
expression $(\ref{mudep})$ to obtain the predicted values.

  2. (b)

Seminaive algorithm

Here, we follow the method in Schelin and Sjöstedt-de Luna (2014). Let $C_{i}$
denote the indicator variable, indicating the presence ($C_{i}=1$) or absence
($C_{i}=0$) of a censored observation. Then we proceed as follows:

    1. 1.

Set $\hat{\mathbf{Z}}_{0}=(\mathbf{Z}_{obs}^{T},\hat{\mathbf{v}}_{0}^{T})$
where $\hat{\mathbf{v}}_{0i}=0$ if $C_{i}=1$.

    2. 2.

Obtain
$\hat{\mbox{\boldmath$\theta$}}_{0}=(\hat{\mbox{\boldmath$\beta$}},\hat{\sigma^{2}},\hat{\phi})$
from $\hat{\mathbf{Z}}_{0}$ by least squares.

    3. 3.

Set $k=0$. Let $\hat{\mathbf{Z}}^{-i}_{k}$ denote the data vector
$\hat{\mathbf{Z}}_{k}$ where unit $i$ is removed.

    4. 4.

Find the predictor $\hat{\mathbf{Z}}(h,\hat{\mbox{\boldmath$\theta$}_{k}})$
for all $i$ such that $C_{i}=1$ using expression (10) and
$\hat{\mathbf{Z}}^{-i}_{k}$, where $h$ denotes the distance between
coordinates (isotropy).

    5. 5.

Set $\hat{\mathbf{Z}}_{k+1}=(\mathbf{Z}_{obs}^{T},\hat{\mathbf{v}}_{k+1}^{T})$
where
$\hat{\mathbf{v}}_{k+1}=(\max(0,\min(\hat{\mathbf{Z}}(h,\hat{\mbox{\boldmath$\theta$}_{k}}),LOD)):C_{i}=1)$.

    6. 6.

Update $\mbox{\boldmath$\theta$}_{k}$ from the new data
$\hat{\mathbf{Z}}_{k+1}$.

    7. 7.

For iteration $k+1$, repeat steps 4-6 until the convergence criterion is
satisfied for some constants $c_{1}$, $c_{2}$, and $c_{3}$ chosen by the user.

    8. 8.

Once convergence is attained, evaluate the estimate
$\hat{\mbox{\boldmath$\theta$}}$ in expression $(\ref{mudep})$ to obtain the
predicted values.

Note that this algorithm only works for the left-censored case. The constants
$c_{1}$, $c_{2}$, and $c_{3}$ are chosen such that the following three
conditions are satisfied:

$\displaystyle\left|\frac{\hat{\sigma^{2}}(\hat{\mathbf{Z}}_{k+1})-\hat{\sigma^{2}}(\hat{\mathbf{Z}}_{k})}{\hat{\sigma^{2}}(\hat{\mathbf{Z}}_{k})}\right|\leq
c_{1},$ $\displaystyle\hat{\sigma^{2}}(\hat{\mathbf{Z}}_{k+1})\leq
c_{2}\hat{\sigma^{2}}(\mathbf{Z}_{obs})\mbox{\quad and}$
$\displaystyle\hat{\xi}{(\hat{\mathbf{Z}}_{k+1}})>c_{3}\hat{\xi}{(\mathbf{Z}_{obs})}$
(11)

where $\xi(X)$ represents the skewness of $X$.

  3. (c)

SAEM algorithm

Here, we proceed as follows:

    1. 1.

Obtain the estimates of the mean and covariance structure by SAEM procedure.
For further details, see Ordonez et al. (2018).

    2. 2.

Impute the censored observations by the approximate first moments
$\hat{\mathbf{Z}}^{(k)}$, obtained from the SAEM procedure.

    3. 3.

Evaluate the above estimates in expression (10) to obtain predicted values for
the unobserved locations. In this case, we use $\mathbf{Z}_{obs^{*}}$ to
denote the observed data and estimate the censored observations via SAEM. This
way, we can easily differentiate it from the fully observed response
$\mathbf{Z}_{obs}$.

###  3.3 Application: Estimation and prediction tools

Now, we exemplify the estimation and prediction methods described previously
and implemented in the CensSpatial package. We use the rainfall dataset from
Parana state (Brazil) analyzed in Diggle and Ribeiro Jr (2002). We carried out
some performance comparison for the four prediction methods by using a cross-
validation study with the first $100$ observations for estimation and the
remaining $43$ observations for comparison of the predictive power of the
estimated methods.

The rainfall dataset does not contain censored data, so left censoring was
generated artificially for the first $100$ observations, following Schelin and
Sjöstedt-de Luna (2014), i.e., we fixed the censoring level at $\alpha$,
sorting the data (in ascending order) and setting the LOD to the 100$\alpha\%$
percentile of the generated sample, i.e., the smallest $\alpha n$ value of
$Z(\mathbf{s}_{1}),\ldots,Z(\mathbf{s}_{n})$. This dataset, called
paranacens25, is available in the package CensSpatial with 25$\%$ censoring
level. The following instructions call the data, then observations that will
be used for estimation are stored in the object dataest and those that will be
used for prediction are stored in the object datapred.

    > data(paranacens25)
    > dataest=paranacens25[1:100,]
    > datapred=paranacens25[101:143,]

In order to specify the mean and covariance components for the SCL model, we
use some descriptive graphs. Figure 1 depicts the rainfall dataset, where the
left panel shows the average rainfall in each coordinate and the right panel
shows the sample variogram for the observed data considering a Gaussian
covariance structure. From this figure, we can conclude that it is plausible
to assume a linear trend of the rainfall with respect to its coordinates. The
behavior showed in the sample variogram seems to indicate that it is also
plausible to consider a Gaussian covariance structure for the stochastic
component. As an additional exercise, we compare Matérn covariance structures
with different $\kappa$ and the Gaussian structure. All these fits present
similar results, but we choose the second one because it is simpler and has an
easier interpretation.

Figure 1: Rainfall dataset: (a) Average rainfall over different locations (b)
Sample variogram and Gaussian correlation structure (solid line).

We fitted the left SCL model defined in (1)-(4) with the specifications
previously discussed in the descriptive analysis. To carry out the fit, we
used the object dataest defined previously and the four estimation methods
available in the CensSpatial package: Seminaive, Naive 1 Naive 2 and SAEM
algorithms. For the Seminaive algorithm, we chose the tuning parameter as
$a_{1}=0.1$, $a_{2}=2$ and $a_{3}=5$ for the convergence criterion as
suggested by Schelin and Sjöstedt-de Luna (2014). For the SAEM procedure, we
performed a Monte Carlo simulation with a sample of size $M=15$, a maximum
number of iterations $W=200$ and a cutoff point $c=0.2$. We used the functions
variofit and likfit from the geoR package to define the initial values for the
SAEM algorithm. These values were used for the remaining algorithms in order
to maintain the same conditions in the comparison of the four methods.

The estimation process via SAEM can be performed using the SAEMSCL function.
Now, we detail how to call the routine and show a print summary of the results
obtained through the est object. Arguments of the function are organized in
such a way that the first line contains the parameters related to the censored
observations, the second one the covariance structure, the third and fourth
ones the model specification and finally the last one refers to the SAEM
settings. The code below also provides a print summary with the parameter
estimates, the linear trend, the covariance structure considered, some
information criteria and specific details about the SAEM procedure. Finally,
the output for this function is a list with all the information contained in
the print summary and specific information of the estimation procedure as the
estimates obtained in each iteration. The SAEM parameters for this function
must not be changed unless the user knows how they work. We refer to Lachos et
al. (2017) for details.

    > cov.ini=c(900,68); y=dataest$y1; cc=dataest$cc
    > coords=dataest[,1:2]; coordspred=datapred[,1:2]

    > est=SAEMSCL(cc,y,cens.type="left",trend="1st",coords=coords,
    + kappa=2, cov.model="gaussian", fix.nugget=F,nugget=170,
    + inits.sigmae=cov.ini[2],inits.phi=cov.ini[1], search=T,
    + lower=c(0.0001,0.0001),upper=c(50000,5000),
    + M=15,perc=0.25,MaxIter=200,pc=0.2)

    -------------------------------------------------------------------------
      Spatial Censored Linear regression with Normal errors (SAEM estimation)
    -------------------------------------------------------------------------

    *Type of trend: Linear function of its coordinates,
    mu = beta0 + beta1*CoordX + beta2*CoordY

    *Covariance structure: Gaussian
    ---------
    Estimates
    ---------
           Estimated
    beta 0  367.1364
    beta 1   -0.0758
    beta 2   -0.3124
    sigma2  832.6448
    phi     211.2341
    tau2    360.1578
    ------------------------
    Model selection criteria
    ------------------------
            Loglik     AIC     BIC AICcorr
    Value -353.601 719.202 734.833 720.105
    -------
    Details
    -------
    Type of censoring = left
    Convergence reached? = FALSE
    Iterations = 200 / 200
    MC sample = 15
    Cut point = 0.2

For the Seminaive, Naive 1 and Naive 2 procedures, the estimation process can
be carried out using the functions Seminaive and algnaive12 ( Naive 1 and
Naive 2 algorithms) respectively. As for the SAEMSCL routine, the arguments
are organized in such a way that parameters data, y.col, coords.col, covar, cc
and trend are related to the censored data structure; parameters cov.model,
fix.nugget, nugget and kappa refer to the covariance structure and thetaini,
cons and MaxIter refer to the algorithm specifications as initial values and
stopping criterion. The print summary and the output for these functions are
quite similar to those offered by the SAEMSCL routine. For our example, the
output generated by this function was saved in the objects r and s for
Seminaive and Naive 1, Naive 2 algorithms, respectively, as can be seen in the
following instructions:

    > r=Seminaive(data=dataest,y.col=3,coords.col=1:2,covar=F,cc=cc,trend="1st",
    + cov.model="gaussian",fix.nugget=F,nugget=170,kappa=2,
    + thetaini=c(900,68),cons=c(0.1,2,0.5),MaxIter=200,copred=coordspred)

    > s=algnaive12(data=dataest,y.col=3,coords.col=1:2,covar=F,cc=cc,trend="1st",
    + cov.model="gaussian",fix.nugget=F,nugget=170,kappa=2,
    + thetaini=c(900,68),copred=coordspred)

Table 1 shows the estimation results obtained through the three methods. The
SAEMSCL method (in boldface) produced the best fit in terms of AIC and BIC
values. The second best fitted method was Naive 1.

Table 1: Rainfall dataset. Models fitted for each spatial estimation and
prediction method

Method | $\widehat{\beta}_{0}$ | $\widehat{\beta}_{1}$ | $\widehat{\beta}_{2}$ | $\widehat{\sigma}^{2}$ | $\widehat{\phi}$ | $\widehat{\tau}^{2}$ | $Loglik$ | $AIC$ | $BIC$
---|---|---|---|---|---|---|---|---|---
Seminaive | 354.8862 | -0.0710 | -0.2261 | 771.4669 | 198.1979 | 277.8797 | -439.38 | 890.76 | 906.391
Naive1 | 354.2134 | -0.0671 | -0.2383 | 878.7536 | 119.5233 | 232.2021 | -438.905 | 889.811 | 905.442
Naive2 | 452.3357 | -0.2091 | -0.5186 | 1969.1575 | 49.1863 | 579.2321 | -501.812 | 1015.625 | 1031.256
SAEM | 367.1364 | -0.0758 | -0.3124 | 832.6448 | 211.2341 | 360.1578 | -353.601 | 719.202 | 734.833

Now we proceed to make predictions under each algorithm to compare their
performance. To obtain predictions via SAEM, we call the predSCL function,
which has three arguments: xpred, which indicates the covariate values at
prediction locations; the prediction locations coordspred; and est, which is
an object of the class "SAEMSpatialCens" and it is returned by the SAEMSCL
function. This routine provides a list with the predicted values along with
their standard deviations. The following instructions save the output of this
function in the object h for the rainfall averages case. The prediction
coordinates were specified using the object coordspred defined previously.

    > xpred=cbind(1,coordspred)
    > h=predSCL(xpred=xpred,coordspred = coordspred,est = est)

The functions Seminaive and algnaive12 already provide the predicted values
for a set of locations copred of interest. For the rainfall average data, if r
and s are objects of type Seminaive and algnaive12, the predicted values are
provided as r$predictions, s$predictions1 and s$predictions2.

We can use the results given by these functions to construct graphs that
facilitate the prediction interpretation, e.g., panel (b) of Figure 2 was
plotted using the output of the predSCL function following the next code
lines:

    > reval=datapred$y1; predsaem= h$prediction; sdsaem=h$sdpred
    > LI=predsaem - (1.96*sdsaem); LS=predsaem + (1.96*sdsaem)

    > plot(reval,type="l",ylim=c(100,450),ylab="Rainfall Average",xlab="")
    > lines(predsaem,col="gray48",type="both",pch=20)
    > lines(LI,col="gray10",lty=2)
    > lines(LS,col="gray10",lty=2)
    > legend(locator(1),c("Real Data", "Prediction SAEM", "Confidence bands"),
    + col=c(1,"gray48","gray10"),lty=c(1,1,2),pch=c(NA,20,NA),bty="n")

Figure 2 shows the predicted values at the observed locations for the four
methods. In the left panel (a) we can compare the methods with the true
observed values in solid lines. All methods seem to follow the data properly.
The right panel (b) shows the real data and a 95$\%$ confidence interval for
the prediction obtained via SAEM.

Figure 2: Rainfall dataset: (a) Spatial predictions at the observed locations
(b) Real data and a 95% confidence interval for the prediction via SAEM
method.

Predictions for the SAEM algorithm may also be assessed graphically using the
predgraphics function. This function has grid1 as argument, which consists of
a grid of coordinates to construct the plot for prediction, xpred, which are
the covariate values defined for this grid, est is an object of the class
"SAEMSpatialCens" and sdgraph is a logical argument which is TRUE by default,
providing the standard deviation of the predictors. The following instructions
illustrate how to use this routine (predgraphics) for the rainfall data.

    > coorgra1=seq(min(coords[,1]),max(coords[,1]),length=50)
    > coorgra2=seq(min(coords[,2]),max(coords[,2]),length=50)

    > grid1=expand.grid(x=coorgra1,y=coorgra2)
    > xpred=cbind(1,grid1)

    > graf=predgraphics(xpred=xpred,est=est,grid1=grid1,points=T,sdgraph=T,
    + colors = gray.colors(300),obspoints = sample(1:sum(cc==0),50),main2="")

Figure 3 is generated by the predgraphics function, i.e., two intensity plots,
one for the predictors and another for the standard deviation. It also has
options for including observed points for comparison purposes.

Figure 3: Rainfall dataset: (a) Observed values and intensity plot for the
predicted values (b) Intensity plot for the standard deviation of prediction
via SAEM method.

We used the mean square prediction error (MSPE) as a measure of prediction
power. This is defined as (see, Fridley and Dixon, 2007):

$\displaystyle MSPE$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}(Z_{i}-\hat{Z}_{i})^{2},$ (12)

where $Z_{i}$ is the observed value, $\hat{Z}_{i}$ is the predicted value and
$n$ is the number of samples predicted. Since the MSPE is an error measure,
lower values lead to better models. Thus, using the ${MSEP}$ we can also
conclude that the SAEM method presents the best fit, with a $\sqrt{MSEP}$
equal to $5.649$, compared with its competitors, all with $\sqrt{MSEP}$
greater than $5.9$.

##  4 Local influence diagnostic utilities

In this section, we describe all measures used for local diagnostics. All
calculations are based in the Hessian matrix of log-likelihood conditional
expectation and the SAEM estimates obtained through the SAEMSCL routine. We
used two simulated processes to illustrate how the package works to detect
influential observations.

Following Zhu and Lee (2001), we compute the Hessian matrix
$\tfrac{\partial^{2}Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\hat{\mbox{\boldmath$\theta$}})}{\partial\mbox{\boldmath$\theta$}\partial\mbox{\boldmath$\theta$}^{\top}}$
in order to obtain local diagnostic measures for a particular perturbation
scheme as described in later subsections. Let
$\mbox{\boldmath$\theta$}=(\mbox{\boldmath$\beta$}^{\top},\mbox{\boldmath$\alpha$}{{}^{\top}})^{\top}$
with $\mbox{\boldmath$\alpha$}=(\alpha_{1},\alpha_{2},\alpha_{3})^{\top}$
where $\alpha_{1}=\sigma^{2}$, $\alpha_{2}=\phi$ and $\alpha_{3}=\tau^{2}$.
Then, the necessary expressions to compute the Hessian matrix
$\ddot{Q}(\mbox{\boldmath$\theta$}|\hat{\mbox{\boldmath$\theta$}})$ are given
by:

$\displaystyle\frac{\partial^{2}Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\hat{\mbox{\boldmath$\theta$}})}{\partial\mbox{\boldmath$\beta$}\partial\mbox{\boldmath$\beta$}^{\top}}$
$\displaystyle=\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}^{-1}\mathbf{X},$ (13)
$\displaystyle\frac{\partial^{2}Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\hat{\mbox{\boldmath$\theta$}})}{\partial\mbox{\boldmath$\beta$}\partial\alpha_{k}}$
$\displaystyle=\mathbf{X}^{\top}\frac{\partial\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}}\widehat{\mathbf{Z}}-\mathbf{X}^{\top}\frac{\partial\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}}\mathbf{X}\mbox{\boldmath$\beta$},$
(14)
$\displaystyle\frac{\partial^{2}Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\hat{\mbox{\boldmath$\theta$}})}{\partial\alpha_{k}\partial\alpha_{l}}$
$\displaystyle=\mathrm{tr}\left(\frac{\partial\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{l}}\frac{\partial\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}}+\frac{\partial^{2}\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}\alpha_{l}}\mbox{\boldmath$\Sigma$}^{-1}\right)$
(15)
$\displaystyle-\frac{1}{2}\left[\mathrm{tr}\left(\widehat{\mathbf{Z}\mathbf{Z}}^{\top}\frac{\partial^{2}\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}\alpha_{l}}\right)-2\mbox{\boldmath$\beta$}^{\top}\mathbf{X}^{\top}\frac{\partial^{2}\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}\alpha_{l}}\widehat{\mathbf{Z}}+\mbox{\boldmath$\beta$}^{\top}\mathbf{X}^{\top}\frac{\partial^{2}\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}\alpha_{l}}\mathbf{X}\mbox{\boldmath$\beta$}\right].$

###  4.1 Local influence

We derive the normal curvature of the local influence for some common
perturbation schemes that could be present either in the model or the data. We
will consider perturbations in response, scale matrix and explanatory
variables.

Consider a perturbation vector $\omega$ varying in an open region
$\mbox{\boldmath$\Omega$}\subset\mathbb{R}^{n}$. The conditional log-
likelihood expectation of the complete data for the perturbed model will be
denoted by
$Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\hat{\mbox{\boldmath$\theta$}})=E(\ell(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\mathbf{Z}_{c})|\mathbf{V},\mathbf{C},\hat{\mbox{\boldmath$\theta$}})$
which reaches its maximum at
$\hat{\mbox{\boldmath$\theta$}}(\mbox{\boldmath$\omega$})$. Assume there is an
$\mbox{\boldmath$\omega$}_{0}$ in $\Omega$ such that
$\ell(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}_{0}|\mathbf{Z}_{c})=\ell(\mbox{\boldmath$\theta$}|\mathbf{Z}_{c})$
for all $\theta$. We use the $Q$-displacement function defined as:

$f_{Q}(\mbox{\boldmath$\omega$})=2[Q(\mbox{\boldmath$\theta$}|\hat{\mbox{\boldmath$\theta$}})-Q(\mbox{\boldmath$\theta$}|\hat{\mbox{\boldmath$\theta$}}(\mbox{\boldmath$\omega$}))],$

to obtain the influence graph
$\delta(\mbox{\boldmath$\omega$})=(\mbox{\boldmath$\omega$}^{\top},f_{Q}(\mbox{\boldmath$\omega$}))^{\top}$.
The local behavior of $f_{Q}(\mbox{\boldmath$\omega$})$ can be summarized
using the curvature $C_{f_{Q,\mathbf{d}}}$ at $\mbox{\boldmath$\omega$}_{0}$
in the direction of some unit vector $\mathbf{d}$. This curvature is defined
by:

$C_{f_{Q,\mathbf{d}}}=-2\mathbf{d}^{\top}\ddot{Q}_{\mbox{\boldmath$\omega$}_{0}}\mathbf{d}\hskip
5.69046pt\quad\textup{and}\quad\hskip
5.69046pt-\ddot{Q}_{\mbox{\boldmath$\omega$}_{0}}=\Delta_{\mbox{\boldmath$\omega$}_{0}}^{\top}-\ddot{Q}(\mbox{\boldmath$\theta$}|\hat{\mbox{\boldmath$\theta$}})^{-1}\Delta_{\mbox{\boldmath$\omega$}_{0}},$

where
$\Delta_{\mbox{\boldmath$\omega$}}=\frac{\partial^{2}Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\hat{\mbox{\boldmath$\theta$}})}{\partial\mbox{\boldmath$\theta$}\partial\mbox{\boldmath$\omega$}^{\top}}$
and $\ddot{Q}(\mbox{\boldmath$\theta$}|\hat{\mbox{\boldmath$\theta$}})$ is the
Hessian matrix with elements as in expressions (13)-(13).

As described in Cook (1986), the quantities $C_{f_{Q,\mathbf{d}}}$ and
$\ddot{Q}_{\mbox{\boldmath$\omega$}_{0}}$ are useful for detecting influential
observations. Let the spectral (orthonormal) decomposition of
$-2\ddot{Q}_{\mbox{\boldmath$\omega$}_{0}}=\sum_{i=1}^{n}{\lambda_{i}\mathbf{a}_{i}\mathbf{a}_{i}^{\top}}$
with
$\lambda_{1}\geq\lambda_{2}\geq\cdots\lambda_{r}>\lambda_{r+1}=\cdots=\lambda_{n}=0$
eigenvalues and $\mathbf{a}_{i},$ for $i=1,\cdots,n$, eigenvectors of this
matrix. Zhu and Lee (2001) proposed inspecting all eigenvectors corresponding
to nonzero eigenvalues to capture the relevant information about influential
observations.

Consider the aggregated contribution of each eigenvector with a corresponding
nonzero eigenvalue as
$\tilde{\lambda_{i}}=\lambda_{i}/(\lambda_{1}+\cdots+\lambda_{r})$. Let
$\mathbf{a}_{i}^{2}=(a_{11}^{2},\cdots,a_{in}^{2})^{\top}$ and
$M(0)=\sum_{i=1}^{r}\lambda_{i}\mathbf{a}_{i}^{2}$ with
$M(0)_{l}=\sum_{i=1}^{r}\lambda_{i}\mathbf{a}_{il}^{2}$ for $l=1,\cdots,n$,
the $l$-th component of $M(0)$. The assessment of influential observations can
be performed through visual inspection of $M(0)_{l}$ plotted against the index
$l$. The $l$-th case may be considered influential if $M(0)_{l}$ is larger
than a benchmark value. Zhu and Lee (2001) also considered the conformal
normal curvature
$B_{f_{Q,\mathbf{d}}}(\mbox{\boldmath$\theta$})=C_{f_{Q,\mathbf{d}}}(\mbox{\boldmath$\theta$})/tr[-2\ddot{Q}_{\mbox{\boldmath$\omega$}_{0}}]$,
whose computation is quite simple and has the property that $0\leq
B_{f_{Q,\mathbf{d}}}(\mbox{\boldmath$\theta$})\leq 1$. They also showed that
for a perturbation vector $\mathbf{d}_{l}$ whose $l$-th entry is 1 and the
remaining entries are 0,
$M(0)_{l}=B_{f_{Q,\mathbf{d}_{l}}}(\mbox{\boldmath$\theta$})$, so we can
obtain $M(0)_{l}$ through
$B_{f_{Q,\mathbf{d}_{l}}}(\mbox{\boldmath$\theta$})$.

In general, there is no rule to judge how influential a specific case in the
data is. Some measures have been proposed to determine the benchmark value in
the visual inspection of $M(0)$. Let $\bar{M}(0)$ and $SM(0)$ be the mean and
standard error of the vector $M(0)$ respectively. Poon and Poon (1999)
proposed using $2\bar{M}(0)$ as a benchmark value while Zhu and Lee (2001)
proposed $\bar{M}(0)+2SM(0)$, for taking into account the standard deviation
of $M(0)$. In the CensSpatial package, we use the benchmark proposed by Lee
and Xu (2004) as $\bar{M}(0)+c^{*}SM(0)$, where $c^{*}$ is a constant, which
can vary depending of the specific application.

###  4.2 Perturbation schemes

We consider three types of perturbation schemes: perturbation in the response
variable, which can indicate influence in their own predicted values,
perturbation in the scale matrix $\Sigma$, which may reveal the individuals
that are influential in the estimation of the scale parameter $\alpha$ and
finally perturbation in explanatory variables. For all these schemes, we are
interested in computing the $\Delta_{\mbox{\boldmath$\omega$}}$ matrix with
elements:

$\Delta_{\mbox{\boldmath$\beta$}}=\frac{\partial^{2}Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\hat{\mbox{\boldmath$\theta$}})}{\partial\mbox{\boldmath$\beta$}\partial\mbox{\boldmath$\omega$}^{\top}}\hskip
2.84544pt\quad\textup{and}\quad\hskip
2.84544pt\Delta_{\mbox{\boldmath$\alpha$}}=(\Delta_{\alpha_{1}}^{\top},\Delta_{\alpha_{2}}^{\top},\Delta_{\alpha_{3}}^{\top}),$

where
$\Delta_{\alpha_{k}}=\frac{\partial^{2}Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\omega$}|\hat{\mbox{\boldmath$\theta$}})}{\partial\alpha_{k}\partial\mbox{\boldmath$\omega$}^{\top}}$,
for $k=1,2,3$.

####  4.2.1 Response perturbation

A perturbation in response variable $\mathbf{V}$ can be introduced by
replacing it with
$\mathbf{V}(\mbox{\boldmath$\omega$})=\mathbf{V}+\mbox{\boldmath$\omega$}$,
recalling the non-perturbed model $V_{i}=V_{0i}$ for $C_{i}=0$ and that
$V_{i}$ belongs to $[V_{1i},V_{2i}]$ for $C_{i}=1$. In this way we can write
the perturbed model as:

$\displaystyle C_{i}=\left\\{\begin{array}[]{ccc}1&if&V_{1i}(\omega_{i})\leq
Z_{i}(\omega_{i})\leq V_{2i}(\omega_{i})\,,\\\
0&if&Z_{i}(\omega_{i})=V_{0i}(\omega_{i})\,,\end{array}\right.$ (18)

where
$\mathbf{V}_{1}(\mbox{\boldmath$\omega$})=\mathbf{V}_{1}+\mbox{\boldmath$\omega$}$,
$\mathbf{V}_{2}(\mbox{\boldmath$\omega$})=\mathbf{V}_{2}+\mbox{\boldmath$\omega$}$
and
$\mathbf{Z}(\mbox{\boldmath$\omega$})=\mathbf{Z}-\mbox{\boldmath$\omega$}$. To
obtain the perturbed function
$Q(\mbox{\boldmath$\theta$}|\hat{\mbox{\boldmath$\theta$}},\mbox{\boldmath$\omega$})$,
we replace $\widehat{\mathbf{Z}}$ and $\widehat{\mathbf{Z}\mathbf{Z}^{\top}}$
with
$\widehat{\mathbf{Z}}_{\omega}=\widehat{\mathbf{Z}}-\mbox{\boldmath$\omega$}$
and
$\widehat{\mathbf{Z}}_{\omega}\widehat{\mathbf{Z}_{\omega}}^{\top}=\widehat{\mathbf{Z}\mathbf{Z}^{\top}}-\widehat{\mathbf{Z}}\mbox{\boldmath$\omega$}-\mbox{\boldmath$\omega$}\widehat{\mathbf{Z}}^{\top}+\mbox{\boldmath$\omega$}\mbox{\boldmath$\omega$}^{\top}$
in the non-perturbed
$Q(\mbox{\boldmath$\theta$}|\hat{\mbox{\boldmath$\theta$}})$ from Subsection
2.3. Let $\mbox{\boldmath$\omega$}_{0}=\mathbf{0}$ be the vector that
represents the case of no perturbation in response. Then, the matrix
$\Delta_{\mbox{\boldmath$\omega$}_{0}}$ has the elements
$\Delta_{\mbox{\boldmath$\beta$}}=-\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}^{-1}$
and $\Delta_{\mbox{\boldmath$\alpha$}}=[\Delta_{\alpha_{ki}}],$ where

$\Delta_{\alpha_{ki}}=\frac{1}{2}\left((\mathbf{Z}_{(i)}+\mathbf{Z}_{(i)}^{\top})\frac{\partial\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}}\right)+\mathbf{Z}_{(i)}^{\top}\mathbf{P}^{\top}\frac{\partial\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}},$

with
$\mathbf{P}=\mathbf{X}(\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}^{-1}\mathbf{X})^{-1}\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}^{-1}$
and $\mathbf{Z}_{(i)}$ the matrix with the $i$-th row being
$\widehat{\mathbf{Z}}$ and the remaining rows being zero. Expressions for
$\tfrac{\partial\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}}$ are
defined in the Appendix.

####  4.2.2 Scale matrix perturbation

Now, we consider the scheme perturbation of the form
$\mbox{\boldmath$\Sigma$}(\mbox{\boldmath$\omega$})=\mathbf{D}(\mbox{\boldmath$\omega$})\mbox{\boldmath$\Sigma$}$
to study the effects of perturbation on the scale matrix. In this case
$\mathbf{D}(\mbox{\boldmath$\omega$})$ is an $n\times n$ diagonal matrix with
elements $\omega$. The non-perturbed model is obtained when $\omega_{i}=1$ for
all $i=1,\cdots,n$ and therefore $\Delta_{\mbox{\boldmath$\omega$}_{0}}$ is a
$(p+3)\times n$ matrix with components
$\Delta_{\mbox{\boldmath$\beta$}_{i}}=\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}^{-1}\mathbf{d}(i)[\mathbf{I_{n}}-\mathbf{P}]\mathbf{Z}$
and $\Delta_{\mbox{\boldmath$\alpha$}}=[\Delta_{\alpha_{ki}}]$, with

$\Delta_{\alpha_{ki}}=-\mathrm{tr}\left(\mbox{\boldmath$\Sigma$}^{-1}\mathbf{d}(i)\frac{\partial\mbox{\boldmath$\Sigma$}}{\partial\alpha_{k}}\right)\\\
-\frac{1}{2}\left[\mathrm{tr}\left(\widehat{\mathbf{Z}\mathbf{Z}^{\top}}\left(\mbox{\boldmath$\Sigma$}^{-1}\frac{\partial\mbox{\boldmath$\Sigma$}}{\partial\alpha_{k}}\mbox{\boldmath$\Sigma$}^{-1}\right)\mathbf{d}(i)\right)+\widehat{\mathbf{Z}}^{\top}[\mathbf{P}^{\top}-2\mathbf{I}]\left(\mbox{\boldmath$\Sigma$}^{-1}\frac{\partial\mbox{\boldmath$\Sigma$}}{\partial\alpha_{k}}\mbox{\boldmath$\Sigma$}^{-1}\right)\mathbf{d}(i)\mathbf{P}\widehat{\mathbf{Z}}\right],$

for $k=1,2,3$, where $\mathbf{d}(i)$ is an $n\times n$ matrix with the $i$-th
diagonal element equal to 1 and the remaining elements equal to zero.

####  4.2.3 Explanatory variable perturbation

We consider for this scheme a perturbation of the form
$\mathbf{X}(\mbox{\boldmath$\omega$})=\mathbf{X}+\mathbf{W}$ and replace the
perturbed
$Q(\mbox{\boldmath$\theta$}|\hat{\mbox{\boldmath$\theta$}},\mbox{\boldmath$\omega$})$
function as in the response perturbation, in this case,
$\mathbf{W}=\mbox{\boldmath$\omega$}\mathbf{1}^{\top}$. Under the non-
perturbed model, $\mbox{\boldmath$\omega$}=\mathbf{0}$ and therefore
$\Delta_{\mbox{\boldmath$\omega$}_{0}}$ has elements
$\Delta_{\mbox{\boldmath$\beta$}_{i}}=\widehat{\mathbf{Z}}^{\top}[\mathbf{I}-\mathbf{2P^{\top}}]\mbox{\boldmath$\Sigma$}^{-1}\mathbf{W_{\textit{{i}}}^{(1)}}$
and $\Delta_{\mbox{\boldmath$\alpha$}}=[\Delta_{\alpha_{ki}}]$, with

$\Delta_{\alpha_{ki}}=\widehat{\mathbf{Z}}^{\top}[\mathbf{I-P^{\top}}]\frac{\partial\mbox{\boldmath$\Sigma$}^{-1}}{\partial\alpha_{k}}\mathbf{W_{\textit{{i}}}^{(1)}}(\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}^{-1}\mathbf{X})^{-1}\mathbf{X}^{\top}\mbox{\boldmath$\Sigma$}\widehat{\mathbf{Z}},$

where $\mathbf{W_{\textit{{i}}}^{(1)}}$ is an $n\times p$ matrix with the
$i$-th row equal to 1 and the remaining rows equal to zero.

###  4.3 Application: Diagnostic tools

In this subsection, we use simulated data to illustrate the local diagnostic
utilities that the CensSpatial package offers. The dataset corresponds to two
linear left censored processes with the same mean structure but with different
covariance structures for errors, i.e., we set a process as in equations
(1)-(4) with spherical and Matérn ($\kappa=0.3$) covariance structures and two
covariates $\mathbf{X}_{1}\sim U(0,1)$ and $\mathbf{X}_{2}\sim U(0,3)$ for the
mean structure. We set a censoring level of $15\%$ and the observations
$\\#91$,$\\#126$ and $\\#162$ as atypical points for the two processes.

To generate the spatial random samples used in the simulations, we call the
function rspacens, its principal arguments are beta and cov.ini to fix the
parameters of the mean and covariance structure respectively, coords to define
the coordinates used for the sample, cens for the censoring level, cens.type
to define the censoring type (left or right) and cov.model to define the
spatial correlation function to use. The code below illustrates how we
generate the samples and the atypical points used in this paper.

    > nobs=200; npred=100
    > r1=sample(seq(1,30,length=400),n+n1); r2=sample(seq(1,30,length=400),n+n1)
    > coords=cbind(r1,r2); cov.ini=c(2,0.1); beta=c(5,3,1)
    > xtot<-cbind(1,runif((n+n1)),runif((n+n1),2,3)); xobs=xtot[1:n,]

    > obj=rspacens(cov.pars=c(3,.3,0),beta=beta,x=xtot,coords=coords,cens=0.15,
    + n=(nobs+ npred),n1=npred,cov.model="matern" ,cens.type="left",kappa=0.3)

    > y=obj$datare[,3]

    > y[91]=y[91]+ 5*sd(y); y[126]=y[126]+ 5*sd(y); y[162]=y[162]+ 5*sd(y)

To assess the influence of outliers in the estimation process, we use the
localinfmeas function. This uses an object of class "SAEMSpatialCens", i.e.,
an object as a result of the SAEM estimation procedure to compute the value of
the vector $M(0)$ for each perturbation scheme and detect the influential
points using the Lee and Xu (2004) criteria. The benchmark used for the
routine is $c^{*}=3$ by default. On the other hand, the function atypical
filters the results to get only the atypical points detected by the
localinfmeas routine. The following instructions generate the output of the
localinfmeas and atypical functions, considering a linear process with Matérn
covariance structure and three perturbation schemes discussed previously.

    > sest=SAEMSCL(cc,y,cens.type="left",trend="other",x=xobs,coords=coords,
    + M=15,perc=0.25,MaxIter=5,pc=0.2,cov.model=type,kappa=0.3,fix.nugget=T,
    + nugget=0,inits.sigmae=cov.ini[1], inits.phi=cov.ini[2],search=T,
    + lower=0.00001,upper=50)

    > w = localinfmeas(sest, fix.nugget = TRUE, c = 3)
    > atypical(w)$RP

                 obs         m0
    32  atypical obs 0.04478987
    91  atypical obs 0.09086417
    126 atypical obs 0.10935255
    162 atypical obs 0.09978107

In these lines, w is a localinfmeas object. Then atypical(w) will be a list
with three elements: EP, SP and RP, containing the atypical detected
observations for each perturbation scheme. By default, localinfmeas also
generates a diagnostic plot with for the three perturbation schemes.

Figure 4: Local influence diagnostic plots for spherical (first row) and
Matérn (second row) processes under (a,d) response perturbation (b,e) scale
matrix perturbation and (c,f) explanatory variable perturbation schemes

Figure 4 (provided by the localinfmeas function) shows the diagnostic plots
for the two simulated datasets. It can be seen that the artificial outliers
were detected, as expected. For the response perturbation, the procedure
detected two atypical observations under the spherical covariance structure
and one under the Matérn covariance structure. Thus, we can conclude that the
proposed local influence measures work quite well to detect atypical
observations in the context of censored spatial models.

##  5 Conclusions

This work introduces an R package for censored spatial data analysis that
allows estimation, prediction and local influence diagnosis for censored
Gaussian processes considering different linear trends and covariance
structures. We discuss in detail all functions available in this package as
well as the theoretical framework. Plots provided by the package constitute
one of the most powerful tools for comparing models, assessing predictions and
analyzing local influence. Applications and several examples can be found in
the manual of the CensSpatial package, available in the CRAN repository. To
the best of our knowledge, this R package is the first proposal, available to
practitioners, to perform statistical analysis of censored spatial data.

Future extensions of this work include the use of scale mixtures of normal
(SMN) distributions to accommodate heavy-tailed features (De Bastiani et al.,
2014). We can also extend the proposal of Diggle et al. (1998) by including
censoring in generalized linear mixed models (GLMMs).

## Acknowledgments

Christian E. Galarza acknowledges support from FAPESP-Brazil (Grant
2015/17110-9).

## References

  * Akaike (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE Transactions on Automatic Control 19, 716–723.
  * Arismendi (2013) Arismendi, J. (2013).  Multivariate truncated moments.  Journal of Multivariate Analysis 117(1), 41–75.
  * Cook (1977) Cook, R. (1977).  Detection of influential observation in linear regression.  Technometrics 19, 15–18.
  * Cook (1986) Cook, R. D. (1986).  Assessment of local influence.  Journal of the Royal Statistical Society, Series B 48, 133–169.
  * De Bastiani et al. (2014) De Bastiani, F., A. H. M. de Aquino Cysneiros, M. A. Uribe-Opazo, and M. Galea (2014).  Influence diagnostics in elliptical spatial linear models.  TEST 24(2), 322–340.
  * Delyon et al. (1999) Delyon, B., M. Lavielle, and E. Moulines (1999).  Convergence of a stochastic approximation version of the EM algorithm.  Annals of Statistics 27(1), 94–128.
  * Dempster et al. (1977) Dempster, A., N. Laird, and D. Rubin (1977).  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Society, Series B 39, 1–38.
  * Diggle and Ribeiro (2007) Diggle, P. and P. Ribeiro (2007).  Model-Based Geostatistics.  Springer Series in Statistics.
  * Diggle and Ribeiro Jr (2002) Diggle, P. and P. Ribeiro Jr (2002).  Bayesian inference in gaussian model-based geostatistics.  Geographical and Environmental Modelling 6, 129–146.
  * Diggle et al. (1998) Diggle, P. J., J. Tawn, and R. Moyeed (1998).  Model-based geostatistics.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 47(3), 299–350.
  * Fridley and Dixon (2007) Fridley, B. L. and P. Dixon (2007).  Data augmentation for a Bayesian spatial model involving censored observations.  Environmetrics 18, 107–123.
  * Genz et al. (2020) Genz, A., F. Bretz, T. Miwa, X. Mi, F. Leisch, F. Scheipl, and T. Hothorn (2020).  mvtnorm: Multivariate Normal and t Distributions.  R package version 1.1-1.
  * Giménez and Galea (2013) Giménez, P. and M. Galea (2013).  Influence measures on corrected score estimators in functional heteroscedastic measurement error models.  Journal of Multivariate Analysis 114, 1–15.
  * Jacqmin-Gadda et al. (2000) Jacqmin-Gadda, H., R. Thiebaut, G. Chene, and D. Commenges (2000).  Analysis of left-censored longitudinal data with application to viral load in HIV infection.  Biostatistics 1, 355–368.
  * Jank (2006) Jank, W. (2006).  Implementing and diagnosing the stochastic approximation EM algorithm.  Journal of Computational and Graphical Statistics 15(4), 803–829.
  * Kuhn and Lavielle (2004) Kuhn, E. and M. Lavielle (2004).  Coupling a stochastic approximation version of EM with an MCMC procedure.  ESAIM: Probability and Statistics 8, 115–131.
  * Lachos et al. (2017) Lachos, V. H., L. A. Matos, T. S. Barbosa, A. M. Garay, and D. K. Dey (2017).  Influence diagnostics in spatial models with censored response.  Environmetrics 28(7), e2464.
  * Lee and Xu (2004) Lee, S. Y. and L. Xu (2004).  R influence analysis of nonlinear mixed-effects models.  Computational Statistics & Data Analysis 45, 321–341.
  * Militino and Ugarte (1999) Militino, A. F. and M. D. Ugarte (1999).  Analyzing censored spatial data.  Mathematical Geology 31(5), 551–561.
  * Ordonez et al. (2018) Ordonez, J. A., D. Bandyopadhyay, V. H. Lachos, and C. R. Cabral (2018).  Geostatistical estimation and prediction for censored responses.  Spatial statistics 23, 109–123.
  * Poon and Poon (1999) Poon, W. and Y. Poon (1999).  Conformal normal curvature and assessment of local influence.  Journal of the Royal Statistical Society, Series B (Statistical Methodology) 61, 51–61.
  * Rathbun (2006) Rathbun, S. L. (2006).  Spatial prediction with left-censored observations.  Journal of Agricultural, Biological, and Environmental Statistics 11(3), 317–336.
  * Schelin and Sjöstedt-de Luna (2014) Schelin, L. and S. Sjöstedt-de Luna (2014).  Spatial prediction in the presence of left-censoring.  Computational Statistics & Data Analysis 74, 125–141.
  * Schwarz (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  Annals of Statistics 6(2), 461–464.
  * Toscas (2010) Toscas, P. J. (2010).  Spatial modelling of left censored water quality data.  Environmetrics 21(6), 632–644.
  * Vaida and Liu (2009) Vaida, F. and L. Liu (2009).  Fast implementation for normal mixed effects models with censored response.  Journal of Computational and Graphical Statistics 18, 797–817.
  * Zhu and Lee (2001) Zhu, H. and S. Lee (2001).  Local influence for incomplete-data models.  Journal of the Royal Statistical Society, Series B (Statistical Methodology) 63, 111–126.
  * Zhu et al. (2001) Zhu, H., S. Lee, B. Wei, and J. Zhou (2001).  Case-deletion measures for models with incomplete data.  Biometrika 88, 727–737.