# Area-level spatio-temporal Poisson mixed models for predicting domain counts
and proportions111Supported by the Instituto Galego de Estatística, by the
grants MTM2017-82724-R and PGC2018-096840-B-I00 of the Spanish Ministerio de
Ciencia, Educación y Universidades and by the Xunta de Galicia (Grupos de
Referencia Competitiva ED431C-2016-015 and Centro Singular de Investigación de
Galicia ED431G/01), all of them through the ERDF.

Miguel Boubeta1, María-José Lombardía2, Francisco Marey-Pérez3, Domingo
Morales4
1 Universidade da Coruña, Spain,
2 Universidade da Coruña, CITIC, Spain,
3Universidad de Santiago de Compostela, Escuela Politécnica Superior de
Ingenierìa, Spain,
4 Universidad Miguel Hernández de Elche, Centro de Investigación Operativa,
Spain.
July 23, 2020

###### Abstract

This paper introduces area-level Poisson mixed models with temporal and SAR(1)
spatially correlated random effects. Small area predictors of the proportions
and counts of a dichotomic variable are derived from the new models and the
corresponding mean squared errors are estimated by parametric bootstrap. The
paper illustrates the introduced methodology with two applications to real
data. The first one deals with data of forest fires in Galicia (Spain) during
2007-2008 and the target is modeling and predicting counts of fires. The
second one treats data from the Spanish living conditions survey of Galicia of
2013 and the target is the estimation of county proportions of women under the
poverty line.
Key words: Forest fires, living conditions survey, small area estimation,
Poisson mixed models, count data, bootstrap, poverty proportion.
AMS subject classification: 62E30, 62J12

##  1 Introduction

This paper introduces statistical methodology for estimating counts, in
particular, problems as diverse as the counts of forest fires and poverty
proportions in areas of Galicia. Galicia is an autonomous community in the
northwest of Spain with an economic activity strongly related to natural
resources. During the last years, a problem that the local government faces
and that directly attacks natural resources are fires. Year after year, the
number of forest fires and the area burned in Galicia has been greater than in
other regions of Spain. The objective is to study the number of wildfires by
forest areas and months, so that the government can know their behavior in
each area and take appropriate measures in each.

Another field of interest to apply this methodology is poverty studies. This
work estimates the poverty map for Galician women, county by county. The
increase in socioeconomic differences between the areas of the inner zone,
which are poorer, and the coastal areas, that have greater development, is
worrying.

Therefore, developing statistical methodologies to model counting data and to
predict domain counts and proportions is important to understand the phenomena
studied (like poverty or arsons) and, consequently, to make decisions about
public policies. This manuscript introduces several extensions of the basic
area-level Poisson mixed model, for better fitting the needs of real data and
giving rise to increasingly complex and realistic models.

When auxiliary variables related to the target count variable are available at
the area level, the Poisson mixed models investigated by Boubeta et al. (2015,
2016, 2017, 2019), the multinomial logit mixed models studied by López-
Vizcaíno et al. (2013, 2015) or the compositional mixed model given by Esteban
et al. (2020) link all the domains to enhance the estimation at a particular
area, that is, they borrow strength from other areas. Their models have random
effects taking into account the between-domain variability that is not
explained by the auxiliary variables, but they assume that the domain random
effects are independent. However, in socioeconomic, environmental and
epidemiological applications, estimates for areas that are spatially close may
be more alike than estimates for areas that are further apart. In fact,
Cressie (1993) shows that not employing spatial models may lead to inefficient
inferences when the auxiliary variables does not explain the spatial
correlation of the study variable.

In small area estimation (SAE), modelling the spatial correlation between data
from different areas allows to borrow even more strength from the areas. This
recommendation was applied to the basic Fay-Herriot model by Singh et al.
(2005). Later, several authors have proposed new spatial area-level linear
mixed models. Petrucci and Salvati (2006), Pratesi and Salvati (2008), Molina
et al. (2009), Marhuenda et al. (2013) and Chandra et al. (2015) consider
linear mixed models (LMM) that extend the Fay-Herriot model. Most of these
papers assume that area effects follow a simultaneously autoregressive process
of order 1 or SAR(1).

In the Bayesian framework, Moura and Migon (2002) and You and Zhou (2011)
consider spatial stationary mixed models, Sugasawa et al. (2015) study an
empirical Bayesian estimation method with spatially non-stationary
hyperparameters for area-level discrete and continuous data having a natural
exponential family distribution. Choi et al. (2011) examine several spatio-
temporal mixed models in small area health data applications and develop new
accuracy measures to assess the recovery of true relative risks. They apply
the spatio temporal models to study chronic obstructive pulmonary disease at
county level in Georgia.

Concerning nonparametric and robust methods, Opsomer et al. (2008) give a
small area estimation procedure using penalized spline regression with
applications to spatially correlated data. Ugarte et al. (2006) and Ugarte et
al. (2010) study the geographical distribution of mortality risk using small
area techniques and penalized splines. Chandra et al. (2012) introduce a
geographical weighted empirical best linear unbiased predictor for a small
area average and give an estimator of its conditional mean squared error
(MSE). Baldermann et al. (2016) describe robust SAE methods under spatial non-
stationarity linear mixed models. Chandra et al. (2017) introduce small area
predictors of counts under a non-stationary spatial model. Chandra et al.
(2018) develop a geographically weighted regression extension of the logistic-
normal and the Poisson-normal generalized linear mixed models (GLMM) allowing
for spatial nonstationarity.

A spatio-temporal extension of the Fay-Herriot model was proposed by Singh et
al. (2005) using the Kalman filtering approach. Under this model, they obtain
a second order approximation to the MSE of the EBLUP. Later, Pereira and
Coelho (2012) and Marhuenda et al. (2013) derive empirical best linear
unbiased predictors under similar models. Specifically, they consider a SAR(1)
spatial correlation structure and an AR(1) process for the temporal component.
Esteban et al. (2012, 2016) present new spatio-temporal model by assuming
AR(1) and MA(1)-correlated random effects. They also propose bootstrap
procedures for estimating the mean squared error and they analyse the
behaviour of the proposed model against other simpler models through several
simulation experiments.

The above cited papers introduce SAE procedures that borrows strength from
spatial or temporal correlations. They mainly apply spatial or temporal
temporal LMMs to the small area estimation setup. However, few of them deals
with empirical best predictors (EBP) under spatial or spatio-temporal GLMMs.
This work partially covers that gap and studies an area-level Poisson mixed
model containing SAR(1) spatially correlated domain effects and independent
time effects. The final target is the estimation of domain counts and
proportions.

This paper is organized as follows. Section 2 introduces the new area-level
Poisson mixed model and several particularizations and limit cases. This
section derives an algorithm to calculate the method of moment estimators of
the model parameters and presents three bootstrap algorithms for testing the
significance of the variance and autocorrelation parameters. Section 3 gives
the empirical best predictor of the Poisson parameter and of the domain and
domain-time random effects. This section proposes a parametric bootstrap
procedure to estimate the MSEs of the EBPs. Section 4 illustrates the
developed methodology in a environmental field. The target is to predict
counts of fires by forest areas of Galicia. Section 5 estimates women poverty
proportions by counties of Galicia. Section 6 collects the main conclusions.
The paper has two appendixes. Appendix A contains the mathematical derivations
of the fitting algorithm. Appendix B investigates the behaviour of the
proposed fitting algorithm and empirically compares the performance of the
plug-in and EBP by means of simulation experiments.

##  2 The model

This section extends the area-level Poisson mixed model, introduced by Boubeta
et al. (2016), to the spatio-temporal context. In particular, the section
introduces a model with SAR(1)-correlated domain random effects and with
independent domain-time random effects. The new model is a generalization of
Model 1, investigated by Boubeta et al. (2017).

Let $D$ and $T$ be the total number of domains and time instants respectively,
with the corresponding indices $d=1,\ldots,D$, and $t=1,\ldots,T$. Consider a
spatio-temporal model with two independent vectors

$\hbox{\boldmath$v$}_{1}=\underset{1\leq d\leq
D}{\mbox{col}}(v_{1,d})\quad\mbox{ and
}\quad\hbox{\boldmath$v$}_{2}=\underset{1\leq d\leq
D}{\mbox{col}}(\underset{1\leq t\leq T}{\mbox{col}}(v_{2,dt})),$

containing the domain and the domain-time random effects respectively. The
model assumes that the vector of domain random effects
$\hbox{\boldmath$v$}_{1}$ is spatially correlated, following a SAR(1) process
with unknown autoregression parameter $\rho$ and known proximity matrix $W$,
i.e.

$\hbox{\boldmath$v$}_{1}=\rho\hbox{\boldmath$W$}\hbox{\boldmath$v$}_{1}+\hbox{\boldmath$u$}_{1},$

where $\hbox{\boldmath$u$}_{1}\sim
N_{D}(\hbox{\boldmath$0$},\hbox{\boldmath$I$}_{D})$, $0$ is the $D\times 1$
zero vector and $\hbox{\boldmath$I$}_{D}$ denotes the $D\times D$ identity
matrix. It also assumes that the matrix
$(\hbox{\boldmath$I$}_{D}-\rho\hbox{\boldmath$W$})$ is non-singular. Then,
$\hbox{\boldmath$v$}_{1}$ can be expressed as

$\hbox{\boldmath$v$}_{1}=(\hbox{\boldmath$I$}_{D}-\rho\hbox{\boldmath$W$})^{-1}\hbox{\boldmath$u$}_{1}.$
(2.1)

For the proximity matrix $W$, we assume that it is row stochastic. Then, the
autoregression parameter $\rho$ is a correlation, $\rho\in(-1,1)$, and is
called spatial autocorrelation parameter. Some of the most used proximity
matrices are based on: (i) common borders, (ii) distances and (iii)
$k$-nearest neighbours. In all cases, the proximity matrix $W$ is obtained
from an original proximity matrix $\hbox{\boldmath$W$}^{0}$ with diagonal
elements equal to zero and remaining entries depending on the employed option.
In option (i), the non diagonal elements of $\hbox{\boldmath$W$}^{0}$ are
equal to 1 when the two domains corresponding to the row and the column
indices are regarded as neighbours and zero otherwise. In Option (ii), the
nondiagonal elements of the proximity matrix $\hbox{\boldmath$W$}^{0}$ are
defined by applying a monotonously decreasing function to the domain
distances; for example, by using the inverse function. Finally, the non
diagonal elements of $\hbox{\boldmath$W$}^{0}$ in option (iii) are 1 if they
correspond to the $k$-nearest neighbours of a given domain and zero otherwise.
For each option, the row standardization is carried out by dividing each entry
of $\hbox{\boldmath$W$}^{0}$ by the sum of the elements in its row.
Consequently, $W$ is row stochastic. Equation (2.1) implies that
$\hbox{\boldmath$v$}_{1}\sim
N_{D}(\hbox{\boldmath$0$},\hbox{\boldmath$\Gamma$}(\rho))$, where

$\hbox{\boldmath$\Gamma$}(\rho)=\big{(}\gamma_{d_{1}d_{2}}(\rho)\big{)}_{d_{1},d_{2}=1,\ldots,D}=\hbox{\boldmath$C$}^{-1}(\rho)$
(2.2)

and
$\hbox{\boldmath$C$}(\rho)=(\hbox{\boldmath$I$}_{D}-\rho\hbox{\boldmath$W$})^{\prime}(\hbox{\boldmath$I$}_{D}-\rho\hbox{\boldmath$W$})$.
Equation (2.1) implies that $\hbox{\boldmath$v$}_{1}=\underset{1\leq d\leq
D}{\mbox{col}}(v_{1,d})\sim
N_{D}(\hbox{\boldmath$0$},\hbox{\boldmath$\Gamma$}(\rho))$, where
$\hbox{\boldmath$\Gamma$}(\rho)$ is given in (2.2). Therefore, the density
function of the domain random effects $\hbox{\boldmath$v$}_{1}$ is

$f_{v}(\hbox{\boldmath$v$}_{1})=(2\pi)^{-D/2}|\hbox{\boldmath$\Gamma$}(\rho)|^{-1/2}\exp\left\\{-\frac{1}{2}\hbox{\boldmath$v$}_{1}^{\prime}\hbox{\boldmath$\Gamma$}^{-1}(\rho)\hbox{\boldmath$v$}_{1}\right\\}.$

Further, it holds that $v_{1,d}\sim N\big{(}0,\gamma_{dd}(\rho)\big{)}$ and
$v_{1,d_{2}}|v_{1,d_{1}}\sim
N\big{(}\mu_{d_{2}|d_{1}},\sigma^{2}_{d_{2}|d_{1}}\big{)}$, where

$\mu_{d_{2}|d_{1}}=\frac{\gamma_{d_{1}d_{2}}(\rho)}{\gamma_{d_{1}d_{1}}(\rho)}\,v_{d_{1}},\quad\sigma^{2}_{d_{2}|d_{1}}=\gamma_{d_{2}d_{2}}(\rho)-\frac{\gamma_{d_{1}d_{2}}^{2}(\rho)}{\gamma_{d_{1}d_{1}}(\rho)}.$

The interaction domain-time random effects, $\hbox{\boldmath$v$}_{2}$, are
assumed to be independent over time, i.e.

$\hbox{\boldmath$v$}_{2d}=\underset{1\leq t\leq T}{\mbox{col}}(v_{2,dt})\sim
N(\hbox{\boldmath$0$},\hbox{\boldmath$I$}_{T}),\quad\hbox{\boldmath$v$}_{2}=\underset{1\leq
d\leq D}{\mbox{col}}(\hbox{\boldmath$v$}_{2d})\sim
N(\hbox{\boldmath$0$},\hbox{\boldmath$I$}_{DT}).$

Then, the join density function of the random effects
$\hbox{\boldmath$v$}_{1}$ and $\hbox{\boldmath$v$}_{2}$ is

$f_{v}(\hbox{\boldmath$v$}_{1},\hbox{\boldmath$v$}_{2})=(2\pi)^{-D(T+1)/2}|\hbox{\boldmath$\Gamma$}(\rho)|^{-1/2}\exp\left\\{-\frac{1}{2}\,\hbox{\boldmath$v$}_{1}^{\prime}\hbox{\boldmath$\Gamma$}^{-1}(\rho)\hbox{\boldmath$v$}_{1}-\frac{1}{2}\,\hbox{\boldmath$v$}_{2}^{\prime}\hbox{\boldmath$v$}_{2}\right\\}.$

The distribution of the target variable $y_{dt}$, conditionally on the random
effects $v_{1,d}$ and $v_{2,dt}$, is

$y_{dt}|v_{1,d},v_{2,dt}\sim\mbox{Poisson}(\mu_{dt}),\quad
d=1,\ldots,D,\,\,t=1,\ldots,T,$

where $\mu_{dt}$ denotes the mean of the Poisson distribution. We assume that
$\mu_{dt}$ can be expressed as $\nu_{dt}p_{dt}$, where $\nu_{dt}$ is a known
natural number. In the application to fire data, we take $\nu_{dt}=1$ and
$\mu_{dt}=p_{dt}$. In the application to poverty data, $y_{dt}$ counts the
number of sampled poor women and $\nu_{dt}$ is the corresponding sample size
of women in domain $d$ and time instant $t$. In that case, $\nu_{dt}$ and
$p_{dt}$ can be interpreted as size and probability parameters respectively.
The advantage of using a Poisson model instead of a binomial model is that we
can avoid the calculation of combinatorial numbers with values outside the
computer range. Besides, by the nature of our poverty data problem, $\nu_{dt}$
takes large values and $p_{dt}$ small values. Then, everything points to a
good behavior of the Poisson model in addition to its computational
advantages.

As $\nu_{dt}$ is assumed to be known, the Poisson parameter, $\mu_{dt}$, is
determined if and only if one knows the parameter $p_{dt}$. In what follows,
we will refer to $p_{dt}$ as target parameter. To define the area-level
Poisson mixed model with SAR(1) spatial domain effects and independent time
effects, we express the natural parameter $\log\mu_{dt}$ in terms of a set of
$p$ covariates, i.e.

$\displaystyle\mbox{Model ST1: }\log\mu_{dt}$ $\displaystyle=\log\nu_{dt}+\log
p_{dt}$
$\displaystyle=\log\nu_{dt}+\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt},\quad
d=1,\ldots,D,\,\,t=1,\ldots,T,$

where $\mu_{dt}=\mathbb{E}[y_{dt}|v_{1,d},v_{2,dt}]$,
$\hbox{\boldmath$x$}_{dt}=\underset{1\leq k\leq
p}{\hbox{col}^{\prime}}(x_{dtk})$ is the row vector of auxiliary variables,
$\hbox{\boldmath$\beta$}=\underset{1\leq k\leq p}{\hbox{col}}(\beta_{k})$ is
the column vector of regression coefficients and $\phi_{1}$ and $\phi_{2}$ are
standard deviation parameters. Further, Model ST1 assumes that the $y_{dt}$’s
are independent conditionally on the random effects $\hbox{\boldmath$v$}_{1}$
and $\hbox{\boldmath$v$}_{2}$. If we define $u_{1,d}=\phi_{1}v_{1,d}$ and
$u_{2,dt}=\phi_{2}v_{2,dt}$, then $\phi_{1}^{2}$ and $\phi_{2}^{2}$ are
variance component parameters for $u_{1,d}$ and $u_{2,dt}$ respectively. They
can be interpreted as variabilities between domains and between time periods
within each domain respectively. Particular or limit cases of Model ST1 are:

  * (i)

Model ST11 if $\phi_{2}=0$, i.e. area-level temporal Poisson mixed model with
SAR(1) spatial domain effects;

  * (ii)

Model T1 if $\rho=0$, i.e. area-level temporal Poisson mixed model with
independent domain effects and independent domain-time effects;

  * (iii)

Model T12 if $\phi_{1}=0$, i.e. area-level temporal Poisson mixed model with
independent domain-time effects;

  * (iv)

Model S1, if $\phi_{2}=0$ and $T=1$, i.e. area-level Poisson mixed model with
SAR(1) spatial domain effects;

  * (v)

Model 1 if $\rho=\phi_{2}=0$ and $T=1$, i.e. area-level Poisson mixed model
with independent domain effects;

  * (vi)

Model 0 if $\rho=\phi_{1}=\phi_{2}=0$ and $T=1$, i.e. area-level Poisson
regression model.

Boubeta et al. (2016, 2017) gives applications of Models 1, T12 and T1 to SAE
problems. This paper presents statistical methodology for Model ST1. The
corresponding procedures (fitting algorithms, predictors or MSE estimators)
for the above cited submodels can be obtained by straightforward
particularizations. More general models, as area-level Poisson mixed model
with SAR(1) spatial domain effects and AR(1) correlated time effects (Model
ST2), are not considered because the number of time instants $T=14$ and $T=1$
are rather small in the applications to fire data and poverty, respectively.
In addition, the first application considers the data of fires for two
consecutive years in the months of April to October. That is the period in
which forest fires occur in Galicia. Therefore, the 14 periods are divided
into two groups of 7 in which an autoregressive or moving average periodic
modeling could be proposed. This poses an excessive complexity to model data
with few temporary instants. For these reasons, we have chosen to only
incorporate the correlation that comes from the auxiliary variables and the
information from the past of the objective and auxiliary variables. That is,
we limit the complexity of the possible models to take into account in Model
ST1.

Under Model ST1, it holds that

$\mathbb{P}(y_{dt}|\hbox{\boldmath$v$})=\mathbb{P}(y_{dt}|v_{dt})=\frac{1}{y_{dt}!}\exp\\{-\nu_{dt}p_{dt}\\}\nu_{dt}^{y_{dt}}p_{dt}^{y_{dt}},$

where
$p_{dt}=\exp\\{\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\\}$
represents the target parameter. The probability function of the response
variable $\hbox{\boldmath$y$}=\underset{1\leq d\leq
D}{\hbox{col}}(\underset{1\leq t\leq T}{\hbox{col}}(y_{dt}))$, conditionally
on the random effects
$\hbox{\boldmath$v$}=(\hbox{\boldmath$v$}_{1},\hbox{\boldmath$v$}_{2})$, is

$\mathbb{P}(\hbox{\boldmath$y$}|\hbox{\boldmath$v$})=\prod_{d=1}^{D}\prod_{t=1}^{T}\mathbb{P}(y_{dt}|\hbox{\boldmath$v$}).$

The marginal probability function of $y$ is

$\mathbb{P}(\hbox{\boldmath$y$})=\int_{\mathbb{R}^{D(T+1)}}\mathbb{P}(\hbox{\boldmath$y$}|\hbox{\boldmath$v$})f_{v}(\hbox{\boldmath$v$}_{1},\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}=\int_{\mathbb{R}^{D(T+1)}}\psi(\hbox{\boldmath$y$},\hbox{\boldmath$v$})\,d\hbox{\boldmath$v$},$

where

$\displaystyle\psi(\hbox{\boldmath$y$},\hbox{\boldmath$v$})$
$\displaystyle=f_{v}(\hbox{\boldmath$v$}_{1},\hbox{\boldmath$v$}_{2})\prod_{d=1}^{D}\prod_{t=1}^{T}\frac{\exp\\{-\nu_{dt}p_{dt}\\}\nu_{dt}^{y_{dt}}p_{dt}^{y_{dt}}}{y_{dt}!}$
$\displaystyle=c(\hbox{\boldmath$y$})|\hbox{\boldmath$\Gamma$}(\rho)|^{-1/2}\exp\left\\{-\frac{1}{2}\,\hbox{\boldmath$v$}_{1}^{\prime}\hbox{\boldmath$\Gamma$}^{-1}(\rho)\hbox{\boldmath$v$}_{1}-\frac{1}{2}\,\hbox{\boldmath$v$}_{2}^{\prime}\hbox{\boldmath$v$}_{2}\right\\}$
$\displaystyle\cdot\exp\left\\{-\sum_{d=1}^{D}\sum_{t=1}^{T}\nu_{dt}\exp\\{\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\\}\right\\}$
$\displaystyle\cdot\exp\left\\{\sum_{k=1}^{p}\big{(}\sum_{d=1}^{D}\sum_{t=1}^{T}y_{dt}x_{dtk}\big{)}\beta_{k}+\phi_{1}\sum_{d=1}^{D}y_{d.}v_{1,d}+\phi_{2}\sum_{d=1}^{D}\sum_{t=1}^{T}y_{dt}v_{2,dt}\right\\},$

with
$c(\hbox{\boldmath$y$})=(2\pi)^{-\frac{D(T+1)}{2}}\prod_{d=1}^{D}\prod_{t=1}^{T}\big{(}\nu_{dt}^{y_{dt}}/y_{dt}!\big{)}$
and $y_{d.}=\sum_{t=1}^{T}y_{dt}$.

The method of moments (MM) can be used to estimate the vector of parameters,
$\hbox{\boldmath$\theta$}=(\hbox{\boldmath$\beta$}^{\prime},\phi_{1},\phi_{2},\rho)^{\prime}$,
of Model ST1. A set of natural equations for applying the MM algorithm is

$\displaystyle\begin{split}0&=f_{k}(\hbox{\boldmath$\theta$})=\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}]x_{dtk}-\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}y_{dt}x_{dtk},\quad
k=1,\ldots,p,\\\
0&=f_{p+1}(\hbox{\boldmath$\theta$})=\frac{1}{D}\sum_{d=1}^{D}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]-\frac{1}{D}\sum_{d=1}^{D}y_{d.}^{2},\\\
0&=f_{p+2}(\hbox{\boldmath$\theta$})=\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}]-\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}y_{dt}^{2},\\\
0&=f_{p+3}(\hbox{\boldmath$\theta$})=\frac{1}{D(D-1)}\sum_{d_{1}\neq
d_{2}}^{D}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d_{1}.}y_{d_{2}.}]-\frac{1}{D(D-1)}\sum_{d_{1}\neq
d_{2}}^{D}y_{d_{1}.}y_{d_{2}.},\end{split}$ (2.3)

The MM estimator of $\theta$ is obtained by solving the system (2.3) of
nonlinear equations. The updating formula of the Newton-Raphson algorithm is

$\hbox{\boldmath$\theta$}^{(m+1)}=\hbox{\boldmath$\theta$}^{(m)}-\hbox{\boldmath$H$}^{-1}(\hbox{\boldmath$\theta$}^{(m)})\hbox{\boldmath$f$}(\hbox{\boldmath$\theta$}^{(m)}),$
(2.4)

where

$\hbox{\boldmath$\theta$}=\underset{1\leq k\leq
p+3}{\hbox{col}}(\theta_{k}),\quad\hbox{\boldmath$f$}(\hbox{\boldmath$\theta$})=\underset{1\leq
k\leq
p+3}{\hbox{col}}(f_{k}(\hbox{\boldmath$\theta$})),\quad\hbox{\boldmath$H$}(\hbox{\boldmath$\theta$})=\left(\frac{\partial
f_{k}(\hbox{\boldmath$\theta$})}{\partial\theta_{\ell}}\right)_{k,\ell=1,\ldots,p+3}.$
(2.5)

Appendix A contains the calculations of the expectations appearing in
$\hbox{\boldmath$f$}(\hbox{\boldmath$\theta$})$ and
$\hbox{\boldmath$H$}(\hbox{\boldmath$\theta$})$. The MM algorithm under Model
ST1 keeps the steps of Algorithm (2.4), replacing $\theta$, $H$ and $f$ for
those given in (2.5). As algorithm seeds for $\beta$, $\phi_{1}$ and
$\phi_{2}$, the algorithm may take the ML estimators under the model with no
spatial correlation ($\rho=0$), i.e. under Model T1. They can be obtained by
using, for example, the function glmer of the R library lme4. For $\rho$, it
may take the Moran’s I measure of spatial autocorrelation

$I=\frac{D}{\sum_{d_{1}=1}^{D}\sum_{d_{2}=1}^{D}w_{d_{1}d_{2}}}\frac{\sum_{d_{1}=1}^{D}\sum_{d_{2}=1}^{D}w_{d_{1}d_{2}}(\tilde{v}_{1,d_{1}}-\tilde{v}_{1})(\tilde{v}_{1,d_{2}}-\tilde{v}_{1})}{\sum_{d=1}^{D}(\tilde{v}_{1,d}-\tilde{v})^{2}},$
(2.6)

where $\tilde{v}_{1,d}$, $d=1,\ldots,D$, are the predicted random effects
under Model T1, $\tilde{v}=\frac{1}{D}\sum_{d=1}^{D}\tilde{v}_{1,d}$ and the
$w_{d_{1}d_{2}}$’s are the elements of the proximity matrix $W$. The
asymptotic variance of the MM estimator under Model ST1 can be approximated by
a similar bootstrap algorithm to that described by Boubeta et al. (2017).

This section ends by presenting three bootstrap algorithms for testing the
significance of the variance parameters, $\phi_{1}$ and $\phi_{2}$, and of the
autocorrelation parameter $\rho$. Algorithm 1 gives a bootstrap procedure to
test the hypothesis $H_{0}:\,\phi_{1}=0$. We test Model T12 against Model ST1.

1:Fit the Model ST1 to data and calculate $\hat{\hbox{\boldmath$\beta$}}$,
$\hat{\phi}_{1}$, $\hat{\phi}_{2}$ and $\hat{\rho}$.

2:Fit the Model T$1_{2}$ to data and calculate
$\hat{\hbox{\boldmath$\beta$}}^{0}$ and $\hat{\phi}_{2}^{0}$.

3:For $b=1,\ldots,B$, do

  * i)

Generate a bootstrap resample under $H_{0}:\phi_{1}=0$, i.e.

$\displaystyle v_{2,dt}^{*(b)}$ $\displaystyle\sim
N(0,1),\,\,p_{dt}^{*(b)}=\exp\\{\hbox{\boldmath$x$}_{dt}\hat{\hbox{\boldmath$\beta$}}^{0}+\hat{\phi}_{2}^{0}v_{2,dt}^{*(b)}\\},$
$\displaystyle y_{dt}^{*(b)}$
$\displaystyle\sim\mbox{Poiss}(\nu_{dt}p_{dt}^{*(b)}),\,\,d=1,\ldots,D,\,t=1,\ldots,T.$

  * ii)

Fit the Model ST1 to the bootstrap data
$(y_{dt}^{*(b)},\hbox{\boldmath$x$}_{dt}),\,d=1,\ldots,D,\,t=1,\ldots,T$, and
calculate $\hat{\hbox{\boldmath$\beta$}}^{*(b)}$, $\hat{\phi}_{1}^{*(b)}$,
$\hat{\phi}_{2}^{*(b)}$ and $\hat{\rho}^{*(b)}$.

4:Calculate the $p$-value:
$p=B^{-1}\\#\\{\hat{\phi}_{1}^{*(b)}>\hat{\phi}_{1}\\}$.

Algorithm 1 A bootstrap test for $H_{0}:\phi_{1}=0$

If the null hypothesis $H_{0}:\phi_{1}=0$ is rejected, the significance of the
autocorrelation parameter can be tested. Algorithm 2 presents a bootstrap
procedure for testing $H_{0}:\rho=0$. We test Model T1 against Model ST1.

1:Fit the Model ST1 to data and calculate $\hat{\hbox{\boldmath$\beta$}}$,
$\hat{\phi}_{1}$, $\hat{\phi}_{2}$ and $\hat{\rho}$.

2:Fit the Model T1 to data and calculate $\hat{\hbox{\boldmath$\beta$}}^{0}$,
$\hat{\phi}_{1}^{0}$ and $\hat{\phi}_{2}^{0}$.

3:For $b=1,\ldots,B$, do

  * i)

Generate a bootstrap resample under $H_{0}:\rho=0$, i.e.

$\displaystyle v_{1,d}^{*(b)}$ $\displaystyle\sim N(0,1),\,v_{2,dt}^{*(b)}\sim
N(0,1),\,\,p_{dt}^{*(b)}=\exp\\{\hbox{\boldmath$x$}_{dt}\hat{\hbox{\boldmath$\beta$}}^{0}+\hat{\phi}_{1}^{0}v_{1,d}^{*(b)}+\hat{\phi}_{2}^{0}v_{2,dt}^{*(b)}\\},$
$\displaystyle y_{dt}^{*(b)}$
$\displaystyle\sim\mbox{Poiss}(\nu_{dt}p_{dt}^{*(b)}),\,\,d=1,\ldots,D,\,t=1,\ldots,T.$

  * ii)

Fit the Model ST1 to the bootstrap data
$(y_{dt}^{*(b)},\hbox{\boldmath$x$}_{dt}),\,d=1,\ldots,D,\,t=1,\ldots,T$, and
calculate $\hat{\hbox{\boldmath$\beta$}}^{*(b)}$, $\hat{\phi}_{1}^{*(b)}$,
$\hat{\phi}_{2}^{*(b)}$ and $\hat{\rho}^{*(b)}$.

4:Calculate the $p$-value:
$p=B^{-1}\\#\\{\lvert\hat{\rho}^{*(b)}\rvert>\lvert\hat{\rho}\rvert\\}$.

Algorithm 2 A bootstrap test for $H_{0}:\rho=0$

Finally, Algorithm 3 gives a bootstrap procedure for testing the null
hypothesis $H_{0}:\phi_{2}=0$. We test Model ST11 against Model ST1.

1:Fit the Model ST1 to data and calculate $\hat{\hbox{\boldmath$\beta$}}$,
$\hat{\phi}_{1}$, $\hat{\phi}_{2}$ and $\hat{\rho}$.

2:Fit the Model ST11 to data and calculate
$\hat{\hbox{\boldmath$\beta$}}^{0}$, $\hat{\phi}_{1}^{0}$ and
$\hat{\rho}^{0}$.

3:For $b=1,\ldots,B$, do

  * i)

Generate a bootstrap resample under $H_{0}:\phi_{2}=0$, i.e.

$\displaystyle\hbox{\boldmath$v$}_{1}^{*(b)}$ $\displaystyle\sim
N_{D}(\hbox{\boldmath$0$},\hbox{\boldmath$\Gamma$}(\hat{\rho}^{0})),\,\,p_{dt}^{*(b)}=\exp\\{\hbox{\boldmath$x$}_{dt}\hat{\hbox{\boldmath$\beta$}}^{0}+\hat{\phi}_{1}^{0}v_{1,d}^{*(b)}\\},$
$\displaystyle y_{dt}^{*(b)}$
$\displaystyle\sim\mbox{Poiss}(\nu_{dt}p_{dt}^{*(b)}),\,\,d=1,\ldots,D,\,t=1,\ldots,T.$

  * ii)

Fit the Model ST1 to the bootstrap data
$(y_{dt}^{*(b)},\hbox{\boldmath$x$}_{dt}),\,d=1,\ldots,D,\,t=1,\ldots,T$, and
calculate $\hat{\hbox{\boldmath$\beta$}}^{*(b)}$, $\hat{\phi}_{1}^{*(b)}$,
$\hat{\phi}_{2}^{*(b)}$ and $\hat{\rho}^{*(b)}$.

4:Calculate the $p$-value:
$p=B^{-1}\\#\\{\hat{\phi}_{2}^{*(b)}>\hat{\phi}_{2}\\}$.

Algorithm 3 A bootstrap test for $H_{0}:\phi_{2}=0$

##  3 The predictors

This section gives the EBP and a plug-in predictor of $p_{dt}$ under Model
ST1. The EBP of $p_{dt}$ is obtained from the corresponding best predictor
(BP) by replacing the vector of model parameters $\theta$ by an estimator
$\hat{\hbox{\boldmath$\theta$}}$. As the MM estimators are consistent, they
are employed for calculating the EBPs. To avoid overflow numerical problems in
the calculation of the exact EBP, this section proposes alternative
approximations.

Let $\hbox{\boldmath$y$}_{d}$ be the response vector within the domain $d$,
i.e. $\hbox{\boldmath$y$}_{d}=\underset{1\leq t\leq T}{\mbox{col}}(y_{dt})$.
The conditional distribution of the response variable $y$, given the random
effects $\hbox{\boldmath$v$}_{1}$ and $\hbox{\boldmath$v$}_{2}$, is

$\mathbb{P}(\hbox{\boldmath$y$}|\hbox{\boldmath$v$}_{1},\hbox{\boldmath$v$}_{2})=\prod_{d=1}^{D}\mathbb{P}(\hbox{\boldmath$y$}_{d}|v_{1,d},\hbox{\boldmath$v$}_{2,d}),\quad\mathbb{P}(\hbox{\boldmath$y$}_{d}|v_{1,d},\hbox{\boldmath$v$}_{2,d})=\prod_{t=1}^{T}\mathbb{P}(y_{dt}|v_{1,d},v_{2,dt}),$

where

$\displaystyle\mathbb{P}(y_{dt}|v_{1,d},v_{2,dt})$
$\displaystyle=\frac{1}{y_{dt}!}\exp\\{-\nu_{dt}p_{dt}\\}\nu_{dt}^{y_{dt}}p_{dt}^{y_{dt}}$
$\displaystyle=c_{dt}\exp\left\\{y_{dt}(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt})-\nu_{dt}\exp\\{\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\\}\right\\}.$

###  3.1 The empirical best predictor

The BP of $p_{dt}$ is the unbiased predictor minimizing the MSE. It is
obtained as the conditional expectation
$\hat{p}_{dt}(\hbox{\boldmath$\theta$})=\mathbb{E}_{\hbox{\boldmath$\theta$}}[p_{dt}|\hbox{\boldmath$y$}]$.
It holds that

$\mathbb{E}_{\hbox{\boldmath$\theta$}}[p_{dt}|\hbox{\boldmath$y$}]=\frac{\int_{\mathbb{R}^{D(T+1)}}p_{dt}\mathbb{P}(\hbox{\boldmath$y$}|\hbox{\boldmath$v$}_{1},\hbox{\boldmath$v$}_{2})f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}}{\int_{\mathbb{R}^{D(T+1)}}\mathbb{P}(\hbox{\boldmath$y$}|\hbox{\boldmath$v$}_{1},\hbox{\boldmath$v$}_{2})f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}}=\frac{N_{dt}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})}{B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})},$
(3.1)

where

$\displaystyle N_{dt}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{D(T+1)}}\exp\\{\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\\}(\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\mathbb{P}(y_{\ell\tau}|v_{1,\ell},v_{2,\ell\tau}))f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}$
$\displaystyle=\int_{\mathbb{R}^{D(T+1)}}\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\exp\Big{\\{}(y_{\ell\tau}+\delta_{d\ell}\delta_{t\tau})(\hbox{\boldmath$x$}_{\ell\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,\ell}+\phi_{2}v_{2,d\tau})$
$\displaystyle-\nu_{\ell\tau}\exp\\{\hbox{\boldmath$x$}_{\ell\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,\ell}+\phi_{2}v_{2,\ell\tau}\\}\Big{\\}}f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2},$
$\displaystyle B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{D(T+1)}}(\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\mathbb{P}(y_{\ell\tau}|v_{1,\ell},v_{2,\ell\tau}))f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}$
$\displaystyle=\int_{\mathbb{R}^{D(T+1)}}\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\exp\Big{\\{}y_{\ell\tau}(\hbox{\boldmath$x$}_{\ell\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,\ell}+\phi_{2}v_{2,d\tau})$
$\displaystyle-\nu_{\ell\tau}\exp\\{\hbox{\boldmath$x$}_{\ell\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,\ell}+\phi_{2}v_{2,\ell\tau}\\}\Big{\\}}f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2},$
(3.2)

and $\delta_{d\ell}$ and $\delta_{t\tau}$ are Kronecker deltas, i.e.
$\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ otherwise. The numerator
$N_{dt}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$ can be expressed in
terms of the denominator $B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$ as
$N_{dt}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})=B(\hbox{\boldmath$y$}+\hbox{\boldmath$e$}_{dt},\hbox{\boldmath$\theta$})$,
where $\hbox{\boldmath$e$}_{dt}=\underset{1\leq\ell\leq
D}{\hbox{col}}(\underset{1\leq\tau\leq
T}{\hbox{col}}(\delta_{d\ell}\delta_{t\tau}))$. The EBP of $p_{dt}$ is
$\hat{p}_{dt}=\hat{p}_{dt}(\hat{\hbox{\boldmath$\theta$}})$.

The EBP calculation involves complex integrals in a high-dimensional space.
The integrals are approximated by using an antithetic Monte Carlo algorithm.
The steps are:

  1. 1.

Generate $\hbox{\boldmath$v$}_{1}^{(s_{1})}\sim
N_{D}\big{(}\hbox{\boldmath$0$},\hbox{\boldmath$\Gamma$}(\hat{\rho})\big{)}$,
$v_{2,\ell\tau}^{(s_{2})}$ i.i.d. $N(0,1)$ and calculate
$\hbox{\boldmath$v$}_{1}^{(S_{1}+s_{1})}=-\hbox{\boldmath$v$}_{1}^{(s_{1})}$,
$v_{2,\ell\tau}^{(S_{2}+s_{2})}=-v_{2,\ell\tau}^{(s_{2})}$,
$s_{1}=1,\ldots,S_{1}$, $s_{2}=1,\ldots,S_{2}$, $\ell=1,\ldots,D$,
$\tau=1,\ldots,T$.

  2. 2.

Approximate the EBP of $p_{dt}$ by
$\hat{p}_{dt}(\hat{\hbox{\boldmath$\theta$}})=\hat{N}_{dt}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})/\hat{B}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})$,
where

$\displaystyle\hat{B}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})$
$\displaystyle=\sum_{s_{1}=1}^{2S_{1}}\sum_{s_{2}=1}^{2S_{2}}\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\exp\left\\{y_{\ell\tau}(\hbox{\boldmath$x$}_{\ell\tau}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,\ell}^{(s_{1})}+\hat{\phi}_{2}v_{2,\ell\tau}^{(s_{2})})\right.$
$\displaystyle-\left.\nu_{\ell\tau}\exp\\{\hbox{\boldmath$x$}_{\ell\tau}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,\ell}^{(s_{1})}+\hat{\phi}_{2}v_{2,\ell\tau}^{(s_{2})}\\}\right\\},\quad\hat{N}_{dt}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})=\hat{B}(\hbox{\boldmath$y$}+\hbox{\boldmath$e$}_{dt},\hat{\hbox{\boldmath$\theta$}}).$
(3.3)

As the relationship between the mean and probability parameters is
$\mu_{dt}=\nu_{dt}p_{dt}$, and $\nu_{dt}$ is known (size parameter), the EBP
of $\mu_{dt}$ is obtained as an immediate consequence of the EBP of $p_{dt}$.
That is to say, the EBP of $\mu_{dt}$ is
$\hat{\mu}_{dt}(\hat{\hbox{\boldmath$\theta$}})=\nu_{dt}\hat{p}_{dt}(\hat{\hbox{\boldmath$\theta$}})$.

The EBP calculations in (2) are computationally demanding. For this reason, we
propose an approximation to the BP of $p_{dt}$ (3.1) under Model ST1. Divide
$y$ and
$\hbox{\boldmath$v$}=(\hbox{\boldmath$v$}_{1}^{\prime},\hbox{\boldmath$v$}_{2}^{\prime})^{\prime}$
into two parts
$(\hbox{\boldmath$y$}_{d}^{\prime},\hbox{\boldmath$y$}_{d-}^{\prime})^{\prime}$
and
$(\hbox{\boldmath$v$}_{d}^{\prime},\hbox{\boldmath$v$}_{d-}^{\prime})^{\prime}$,
where $\hbox{\boldmath$y$}_{d}=\underset{1\leq t\leq T}{\hbox{col}}(y_{dt})$
$\hbox{\boldmath$y$}_{d-}=\underset{1\leq i\leq D,\,i\neq
d}{\hbox{col}}(\hbox{\boldmath$y$}_{i})$,
$\hbox{\boldmath$v$}_{d}=(v_{1,d},\hbox{\boldmath$v$}_{2,d}^{\prime})^{\prime}$
and $\hbox{\boldmath$v$}_{d-}=\underset{1\leq i\leq D,\,i\neq
d}{\hbox{col}}(\hbox{\boldmath$v$}_{i})$. The conditional distribution of $y$,
given $v$, is

$\mathbb{P}(\hbox{\boldmath$y$}|\hbox{\boldmath$v$})=\prod_{i=1}^{D}\mathbb{P}(\hbox{\boldmath$y$}_{i}|\hbox{\boldmath$v$}_{i})=\mathbb{P}(\hbox{\boldmath$y$}_{d}|\hbox{\boldmath$v$}_{d})\prod_{i=1,i\neq
d}^{D}\mathbb{P}(\hbox{\boldmath$y$}_{i}|\hbox{\boldmath$v$}_{i})=\mathbb{P}(\hbox{\boldmath$y$}_{d}|\hbox{\boldmath$v$}_{d})\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d-}).$
(3.4)

The p.d.f. of $\hbox{\boldmath$v$}_{d}$ is

$f(\hbox{\boldmath$v$}_{d})=f(v_{1,d})f(\hbox{\boldmath$v$}_{2,d}),$

where $v_{1,d}\sim N(0,\gamma_{dd}(\rho))$ and $\hbox{\boldmath$v$}_{2,d}\sim
N(\hbox{\boldmath$0$},\hbox{\boldmath$I$}_{T})$. The component
$B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$ of the BP of $p_{dt}$ in
(3.1) can be rewritten by using the decomposition of the conditional
probability given in (3.4), i.e.

$\displaystyle
B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})=\int_{\mathbb{R}^{T+1}}\Big{[}\int_{\mathbb{R}^{(D-1)(T+1)}}\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d-})f(\hbox{\boldmath$v$}_{d-}|\hbox{\boldmath$v$}_{d})\,d\hbox{\boldmath$v$}_{d-}\Big{]}\mathbb{P}(\hbox{\boldmath$y$}_{d}|\hbox{\boldmath$v$}_{d})f(\hbox{\boldmath$v$}_{d})\,d\hbox{\boldmath$v$}_{d}.$

As
$\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d-})f(\hbox{\boldmath$v$}_{d-}|\hbox{\boldmath$v$}_{d})=\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d-},\hbox{\boldmath$v$}_{d})f(\hbox{\boldmath$v$}_{d-}|\hbox{\boldmath$v$}_{d})$,
the inner integral is

$\int_{\mathbb{R}^{(D-1)(T+1)}}\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d-},\hbox{\boldmath$v$}_{d})f(\hbox{\boldmath$v$}_{d-}|\hbox{\boldmath$v$}_{d})\,d\hbox{\boldmath$v$}_{d-}=\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d}).$

Therefore, it holds that

$\displaystyle B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{T+1}}\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d})\mathbb{P}(\hbox{\boldmath$y$}_{d}|\hbox{\boldmath$v$}_{d})f(\hbox{\boldmath$v$}_{d})\,d\hbox{\boldmath$v$}_{d}.$

By applying similar developments as for the component
$N_{dt}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$ of (3.1), it holds that

$\displaystyle N_{dt}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{T+1}}\exp\\{\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\\}\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d})\mathbb{P}(\hbox{\boldmath$y$}_{d}|\hbox{\boldmath$v$}_{d})f(\hbox{\boldmath$v$}_{d})\,d\hbox{\boldmath$v$}_{d}.$

Under the assumption

$\mathbb{P}(\hbox{\boldmath$y$}_{d-}|\hbox{\boldmath$v$}_{d})\approx\mathbb{P}(\hbox{\boldmath$y$}_{d-}),$
(3.5)

the BP of $p_{dt}$, $\hat{p}_{dt}(\hbox{\boldmath$\theta$})$, can be
approximated by

$\hat{p}_{dt}^{a}(\hbox{\boldmath$\theta$})=N_{dt}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})/B_{d}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$}),$

where

$\displaystyle N_{dt}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{T+1}}\exp\\{\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\\}\prod_{\tau=1}^{T}\mathbb{P}(y_{d\tau}|v_{1,d},v_{2,d\tau})f(v_{1,d})f(\hbox{\boldmath$v$}_{2,d})\,dv_{1,d}\,d\hbox{\boldmath$v$}_{2,d}$
$\displaystyle=B_{d}^{a}(\hbox{\boldmath$y$}_{d}+\hbox{\boldmath$\epsilon$}_{t},\hbox{\boldmath$\theta$}),\quad\hbox{\boldmath$\epsilon$}_{t}=\underset{1\leq\tau\leq
T}{\hbox{col}}(\delta_{t\tau}),$ $\displaystyle
B_{d}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}}\prod_{\tau=1}^{T}\Big{[}\int_{\mathbb{R}}\exp\Big{\\{}y_{d\tau}(\hbox{\boldmath$x$}_{d\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,d\tau})$
$\displaystyle-\exp\\{\hbox{\boldmath$x$}_{d\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,d\tau}\\}\Big{\\}}f(v_{2,d\tau})dv_{2,d\tau}\Big{]}f(v_{1,d})dv_{1,d}.$
(3.6)

Then, the EBP of $p_{dt}$, $\hat{p}_{dt}(\hat{\hbox{\boldmath$\theta$}})$, can
be approximated as follows.

  1. 1.

Estimate
$\hat{\hbox{\boldmath$\theta$}}=(\hat{\hbox{\boldmath$\beta$}}^{\prime},\hat{\phi}_{1},\hat{\phi}_{2},\hat{\rho})$.

  2. 2.

For $s_{1}=1,\ldots,S_{1}$, $s_{2}=1,\ldots,S_{2}$, $\ell=1,\ldots,D$,
$\tau=1,\ldots,T$, generate $\hbox{\boldmath$v$}_{1}^{(s_{1})}\sim
N_{D}\big{(}\hbox{\boldmath$0$},\hbox{\boldmath$\Gamma$}(\hat{\rho})\big{)}$,
$v_{2,\ell\tau}^{(s_{2})}$ i.i.d. $N(0,1)$ and calculate
$\hbox{\boldmath$v$}_{1}^{(S_{1}+s_{1})}=-\hbox{\boldmath$v$}_{1}^{(s_{1})}$,
$v_{2,\ell\tau}^{(S_{2}+s_{2})}=-v_{2,\ell\tau}^{(s_{2})}$.

  3. 3.

Calculate
$\hat{p}_{dt}^{a}(\hat{\hbox{\boldmath$\theta$}})=\hat{N}_{dt}^{a}(\hbox{\boldmath$y$}_{d},\hat{\hbox{\boldmath$\theta$}})/\hat{B}_{d}^{a}(\hbox{\boldmath$y$}_{d},\hat{\hbox{\boldmath$\theta$}})$,
$d=1,\ldots,D$, $t=1,\ldots,T$, where

$\displaystyle\hat{B}_{d}^{a}(\hbox{\boldmath$y$}_{d},\hat{\hbox{\boldmath$\theta$}})$
$\displaystyle=\sum_{s_{1}=1}^{2S_{1}}\prod_{\tau=1}^{T}\sum_{s_{2}=1}^{2S_{2}}\exp\left\\{y_{d\tau}(\hbox{\boldmath$x$}_{d\tau}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,d}^{(s_{1})}+\hat{\phi}_{2}v_{2,d\tau}^{(s_{2})})\right.$
$\displaystyle-\left.\nu_{d\tau}\exp\\{\hbox{\boldmath$x$}_{d\tau}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,d}^{(s_{1})}+\hat{\phi}_{2}v_{2,d\tau}^{(s_{2})}\\}\right\\},\quad\hat{N}_{dt}^{a}(\hbox{\boldmath$y$}_{d},\hat{\hbox{\boldmath$\theta$}})=B_{d}^{a}(\hbox{\boldmath$y$}_{d}+\hbox{\boldmath$\epsilon$}_{t},\hat{\hbox{\boldmath$\theta$}}).$

The EBP approximation, $\hat{p}_{dt}^{a}(\hat{\hbox{\boldmath$\theta$}})$,
maintains the expression of the EBP of $p_{dt}$ under the area-level Poisson
mixed model with independent time effects studied by Boubeta et al (2017).
However, the domain random effects, $v_{1,d}$, are now generated according to
a SAR(1) process.

###  3.2 The plug-in predictor

The plug-in predictor of $p_{dt}$ is obtained by replacing, in the theoretical
expression of $p_{dt}$, the model parameters by their estimates and the random
effects by their EBPs, i.e.

$\hat{p}_{dt}^{P}=\exp\\{\hbox{\boldmath$x$}_{dt}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}\hat{v}_{1,d}+\hat{\phi}_{2}\hat{v}_{2,dt}\\}.$
(3.7)

As above, the EBPs of $\hbox{\boldmath$v$}_{1}$ and $\hbox{\boldmath$v$}_{2}$
are obtained from the respective BPs. The BP of $v_{1,d}$ is

$\hat{v}_{1,d}(\hbox{\boldmath$\theta$})=\mathbb{E}_{\hbox{\boldmath$\theta$}}[v_{1,d}|\hbox{\boldmath$y$}]=\frac{N_{1,d}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})}{B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})},$

where $B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$ was defined in (3.1)
and

$\displaystyle N_{1,d}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{D(T+1)}}v_{1,d}(\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\mathbb{P}(y_{\ell\tau}|v_{1,\ell},\hbox{\boldmath$v$}_{2,\ell}))f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}$
$\displaystyle=\int_{\mathbb{R}^{D(T+1)}}\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\exp\Big{\\{}y_{\ell\tau}(\hbox{\boldmath$x$}_{\ell\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,\ell}+\phi_{2}v_{2,\ell\tau})$
$\displaystyle-\nu_{\ell\tau}\exp\\{\hbox{\boldmath$x$}_{\ell\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,\ell}+\phi_{2}v_{2,\ell\tau}\\}\Big{\\}}v_{1,d}f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}.$

The EBP of $v_{1,d}$ is
$\hat{v}_{1,d}=\hat{v}_{1,d}(\hat{\hbox{\boldmath$\theta$}})$ and it can be
approximated as follows.

  1. 1.

Generate $\hbox{\boldmath$v$}_{1}^{(s_{1})}\sim
N_{D}\big{(}\hbox{\boldmath$0$},\hbox{\boldmath$\Gamma$}(\hat{\rho})\big{)}$,
$v_{2,\ell\tau}^{(s_{2})}$ i.i.d. $N(0,1)$ and calculate
$\hbox{\boldmath$v$}_{1}^{(S_{1}+s_{1})}=-\hbox{\boldmath$v$}_{1}^{(s_{1})}$,
$v_{2,\ell\tau}^{(S_{2}+s_{2})}=-v_{2,\ell\tau}^{(s_{2})}$,
$s_{1}=1,\ldots,S_{1}$, $s_{2}=1,\ldots,S_{2}$, $\ell=1,\ldots,D$,
$\tau=1,\ldots,T$.

  2. 2.

Calculate
$\hat{v}_{1,d}(\hat{\hbox{\boldmath$\theta$}})=\hat{N}_{1,d}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})/\hat{B}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})$,
where $\hat{B}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})$ was
defined in (2) and

$\displaystyle\hat{N}_{1,d}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})$
$\displaystyle=\sum_{s_{1}=1}^{2S_{1}}\sum_{s_{2}=1}^{2S_{2}}\prod_{\ell=1}^{D}\big{\\{}1+\delta_{\ell
d}(v_{1,\ell}^{(s_{1})}-1)\big{\\}}\prod_{\tau=1}^{T}\exp\left\\{y_{\ell\tau}(\hbox{\boldmath$x$}_{\ell\tau}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,\ell}^{(s_{1})}+\hat{\phi}_{2}v_{2,\ell\tau}^{(s_{2})})\right.$
$\displaystyle-\left.\nu_{\ell\tau}\exp\\{\hbox{\boldmath$x$}_{\ell\tau}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,\ell}^{(s_{1})}+\hat{\phi}_{2}v_{2,\ell\tau}^{(s_{2})}\\}\right\\}.$

The BP of the domain-time random effects $v_{2,dt}$ is

$\hat{v}_{2,dt}(\hbox{\boldmath$\theta$})=\mathbb{E}_{\hbox{\boldmath$\theta$}}[v_{2,dt}|\hbox{\boldmath$y$}]=\frac{N_{2,dt}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})}{B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})},$

where $B(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$ was defined in (3.1)
and

$\displaystyle N_{2,dt}(\hbox{\boldmath$y$},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{D(T+1)}}v_{2,dt}\big{(}\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\mathbb{P}(y_{\ell\tau}|v_{1,\ell},\hbox{\boldmath$v$}_{2,\ell})\Big{)}f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}$
$\displaystyle=\int_{\mathbb{R}^{D(T+1)}}\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\exp\Big{\\{}y_{\ell\tau}(\hbox{\boldmath$x$}_{\ell\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,\ell}+\phi_{2}v_{2,d\tau})$
$\displaystyle-\nu_{\ell\tau}\exp\\{\hbox{\boldmath$x$}_{\ell\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,\ell}+\phi_{2}v_{2,\ell\tau}\\}\Big{\\}}v_{2,dt}f(\hbox{\boldmath$v$}_{1})f(\hbox{\boldmath$v$}_{2})\,d\hbox{\boldmath$v$}_{1}d\hbox{\boldmath$v$}_{2}.$

The EBP of $v_{2,dt}$ is
$\hat{v}_{2,dt}=\hat{v}_{2,dt}(\hat{\hbox{\boldmath$\theta$}})$ and it can be
approximated as follows.

  1. 1.

Generate $\hbox{\boldmath$v$}_{1}^{(s_{1})}\sim
N_{D}\big{(}\hbox{\boldmath$0$},\hbox{\boldmath$\Gamma$}(\hat{\rho})\big{)}$,
$v_{2,\ell\tau}^{(s_{2})}$ i.i.d. $N(0,1)$ and calculate
$\hbox{\boldmath$v$}_{1}^{(S_{1}+s_{1})}=-\hbox{\boldmath$v$}_{1}^{(s_{1})}$,
$v_{2,\ell\tau}^{(S_{2}+s_{2})}=-v_{2,\ell\tau}^{(s_{2})}$,
$s_{1}=1,\ldots,S_{1}$, $s_{2}=1,\ldots,S_{2}$, $\ell=1,\ldots,D$,
$\tau=1,\ldots,T$.

  2. 2.

Calculate
$\hat{v}_{2,dt}(\hat{\hbox{\boldmath$\theta$}})=\hat{N}_{2,dt}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})/\hat{B}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})$,
where $\hat{B}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})$ was
defined in (2) and

$\displaystyle\hat{N}_{2,dt}(\hbox{\boldmath$y$},\hat{\hbox{\boldmath$\theta$}})$
$\displaystyle=\sum_{s_{1}=1}^{2S_{1}}\sum_{s_{2}=1}^{2S_{2}}\prod_{\ell=1}^{D}\prod_{\tau=1}^{T}\big{[}1+\delta_{d\ell}\delta_{t\tau}(v_{2,\ell\tau}^{(s_{2})}-1)\big{]}\exp\left\\{y_{\ell\tau}(\hbox{\boldmath$x$}_{\ell\tau}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,\ell}^{(s_{1})}+\hat{\phi}_{2}v_{2,\ell\tau}^{(s_{2})})\right.$
$\displaystyle-\left.\nu_{\ell\tau}\exp\\{\hbox{\boldmath$x$}_{\ell\tau}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,\ell}^{(s_{1})}+\hat{\phi}_{2}v_{2,\ell\tau}^{(s_{2})}\\}\right\\}.$

The EBPs of $v_{1,d}$ and $v_{2,dt}$ are computationally demanding and
therefore this section proposes approximations, $\hat{v}_{1,d}^{a}$ and
$\hat{v}_{2,dt}^{a}$, to the EBPs $\hat{v}_{1,d}$ and $\hat{v}_{2,dt}$
respectively. Under the assumption (3.5), the BP of the domain random effects,
$\hat{v}_{1,d}(\hbox{\boldmath$\theta$})$, can be approximated by

$\hat{v}_{1,d}^{a}(\hbox{\boldmath$\theta$})=\frac{N_{1,d}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})}{B_{d}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})},$

where $B_{d}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})$ is given
in (3.1) and

$\displaystyle N_{1,d}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\prod_{\tau=1}^{T}\bigg{[}\int_{\mathbb{R}}\exp\Big{\\{}y_{d\tau}(\hbox{\boldmath$x$}_{d\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,d\tau})$
$\displaystyle-$
$\displaystyle\nu_{d\tau}\exp\\{\hbox{\boldmath$x$}_{d\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,d\tau}\\}\Big{\\}}f(v_{2,d\tau})\,dv_{2,d\tau}\bigg{]}v_{1,d}f(v_{1,d})\,dv_{1,d}.$

The difference between the approximated BP,
$\hat{v}_{1,d}^{a}(\hbox{\boldmath$\theta$})$, and the BP based on Model T1 is
that the distribution of the domain random effects is SAR(1) instead of being
i.i.d. $N(0,1)$.

The BP of $v_{2,dt}$, $\hat{v}_{2,dt}(\hbox{\boldmath$\theta$})$, can be
approximated by

$\hat{v}_{2,dt}^{a}(\hbox{\boldmath$\theta$})=\frac{N_{2,dt}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})}{B_{d}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})},$

where

$\displaystyle N_{2,dt}^{a}(\hbox{\boldmath$y$}_{d},\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}}\prod_{\tau=1}^{T}\bigg{[}\int_{\mathbb{R}}(1+\delta_{t\tau}(v_{2,d\tau}-1))\exp\Big{\\{}y_{d\tau}(\hbox{\boldmath$x$}_{d\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,d\tau})$
$\displaystyle-\nu_{d\tau}\exp\\{\hbox{\boldmath$x$}_{d\tau}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,d\tau}\\}\Big{\\}}f(v_{2,d\tau})\,dv_{2,d\tau}\bigg{]}f(v_{1,d})\,dv_{1,d}.$

In this case, the underlying spatial correlation structure slightly
complicates the expression of the approximated BP,
$\hat{v}_{2,dt}^{a}(\hbox{\boldmath$\theta$})$, with respect to that obtained
under Model T1.

To finish this section, we introduce the synthetic predictors

$\hat{p}_{dt}^{syn}=\exp\\{\hbox{\boldmath$x$}_{dt}\hat{\hbox{\boldmath$\beta$}}\\},\quad\hat{\mu}_{dt}^{syn}=\nu_{dt}\hat{p}_{dt}^{syn},$
(3.8)

which do not require the calculation of the EBPs of the random effects. The
synthetic predictors can give a parsimonious solution in those cases in which
the standard deviation parameters are very small and a set of highly
informative auxiliary variables is available.

###  3.3 MSE estimation

As accuracy measure of a predictor (EBP or plug-in) of $p_{dt}$ under Model
ST1, this section considers the MSE. It proposes estimating the MSE of the
predictor of $p_{dt}$ by using a parametric bootstrap algorithm based on the
bootstrap procedure given by González-Manteiga et al. (2008, 2010), since the
analytical estimation is not feasible computationally. The steps are:

  1. 1.

Fit the model to the sample and calculate
$\hat{\hbox{\boldmath$\theta$}}=(\hat{\hbox{\boldmath$\beta$}},\hat{\phi}_{1},\hat{\phi}_{2},\hat{\rho})$
under Model ST1.

  2. 2.

For $d=1,\ldots,D$, $t=1,\ldots,T$, repeat $B$ times, $b=1,\ldots,B$.

    1. (a)

Generate the bootstrap random effects $\hbox{\boldmath$v$}_{1}^{*(b)}\sim
N_{D}\big{(}\hbox{\boldmath$0$},\hbox{\boldmath$\Gamma$}(\hat{\rho})\big{)}$
and $\\{v_{2,dt}^{*(b)}\\}$ i.i.d. $N(0,1)$.

    2. (b)

Calculate the theoretical bootstrap quantity
$p_{dt}^{*(b)}=\text{exp}\\{\hbox{\boldmath$x$}_{dt}\hat{\hbox{\boldmath$\beta$}}+\hat{\phi}_{1}v_{1,d}^{*(b)}+\hat{\phi}_{2}v_{2,dt}^{*(b)}\\}$.

    3. (c)

Generate the response variables
$y_{dt}^{*(b)}\sim\mbox{Poiss}(\nu_{dt}p_{dt}^{*(b)})$.

    4. (d)

Calculate $\hat{\hbox{\boldmath$\theta$}}^{*(b)}$ and the predictor (EBP or
plug-in)
$\hat{p}_{dt}^{*(b)}=\hat{p}_{dt}^{*(b)}(\hat{\hbox{\boldmath$\theta$}}^{*(b)})$.

  3. 3.

Output:

$mse^{*}(\hat{p}_{dt})=\frac{1}{B}\sum_{b=1}^{B}\big{(}\hat{p}_{dt}^{*(b)}-p_{dt}^{*(b)}\big{)}^{2}.$

##  4 Application to fire data

This section presents an application to real data of wildfires in Galicia
during $2007-2008$. The objective of this study is to analyse the target
variable number of wildfires by forest areas and months. The domains are the
forest areas. For each domain, the data set collects the number of wildfires
by month. The section only takes the months between April to October for being
the months with the greatest number of fires. This is the main reason for not
considering models with time correlated random effects. The total number of
domains and time instants are $D=63$ and $T=14$, respectively.

The response variable at domain $d$, $d=1,\ldots,D$, and time $t$,
$t=1,\ldots,T$, $y_{dt}$, is explained by some auxiliary variables through an
area-level spatio-temporal Poisson mixed model with SAR(1)-correlated domain
and independent domain-time effects. As proximity matrix, $W$, we take the
common border option described in Section 2. The auxiliary variables are
acumRain (accumulated rain in $l/m^{2}$), averTemp (Average air temperature in
${}^{\circ}C$) and cadHold (Number of owners of cadastral parcels). The two
first variables are calculated by averaging the measurements of meteorological
stations by months and forests areas. The third variable is annually updated
in the Land Register of the Spanish Ministry of Agriculture, Fisheries and
Food. See Boubeta et al. (2019) for further information about the auxiliary
variables.

Table 4.1 presents the MM estimates of the fixed effects under Model ST1. This
table suggests that $acumRain$ is protective, since it causes a decrease in
the response variable if it increases and the other variables remain fixed. On
the other hand, $averTemp$ and $cadHold$ are directly related to the response
variable since their signs are positive. The three covariates are significant
at the level $\alpha=0.05$.

Table 4.1: MM estimates of regression parameters under Model ST1. Variable | Est. | s.e. |  $z$-value |  $p$-value
---|---|---|---|---
$Intercept$ | 0.7736 | 0.0870 | 8.8910 |  $<$ 0.001
$acumRain$ | -0.5207 | 0.0682 | -7.6318 |  $<$ 0.001
$averTemp$ | 0.1507 | 0.0661 | 2.2820 | 0.0225
$cadHold$ | 0.2931 | 0.0605 | 4.8466 |  $<$ 0.001

The MM estimates of the standard deviations are $\hat{\phi}_{1}=0.331$ and
$\hat{\phi}_{2}=0.696$. Their $95\%$ percentile bootstrap confidence intervals
are $(0.130,0.467)$ and $(0.503,0.829)$ respectively. The estimate of the
spatial autocorrelation is $\hat{\rho}=0.327$. The Algorithms 1 and 3 are used
to test the significance of the parameters $\phi_{1}$ and $\phi_{2}$. The
obtained bootstrap $p$-values are $0.01$ and $0.00$ respectively. The
conclusion is that both standard deviations are significantly different from
$0$. In addition, as the hypothesis $H_{0}:\phi_{1}=0$ is rejected, the
Algorithm 2 is applied to test $H_{0}:\rho=0$. The obtained bootstrap
$p$-value is $0.00$. Then, based on the bootstrap $p$-values, it is
recommended to use a spatio-temporal Poisson mixed model ST1 to analyse
wildfires in Galicia by forest areas and the considered months in the period
$2007-2008$.

Figure 4.1 plots the Pearson residuals of the synthetic estimator,
$\hat{\mu}_{dt}^{syn}$, under Model 0 (left) and of the EBP $\hat{\mu}_{dt}$
under Model ST1 (right). A clear improvement is achieved when one uses a more
complex model, since the Pearson residuals are closer to 0. In this case, the
EBP of $\mu_{dt}$ under Model ST1 is more competitive than the plug-in
predictor under Model $0$, since the empirical MSEs of Pearson residuals are
$0.454$ and $1.247$ respectively.

|
---|---
Figure 4.1: Pearson residuals of the synthetic estimator under Model 0 (left)
and of the EBP of $\mu_{dt}$ under Model ST1 (right).

Figure 4.2 maps the EBPs of numbers of wildfires by forest areas under the
fitted Model ST1. We group the forest areas by number of fires with cuts in 5,
10 and 15. The legends present, between parentheses, the number of forest
areas in each group. We take the same months between August and October. The
figure suggests that the highest number of wildfires are concentrated in
western coastal areas and in the south of the region. On the other hand,
regarding the temporal behaviour, the highest number of fires is found in the
months of $2007$ (specially in September and October, being a year with the
summer delayed), while in $2008$ there was an impressive decrease, because it
was a particularly rainy summer (information contrasted with MeteoGalicia).

---
Figure 4.2: Estimated fires from August to October in 2007-2008.

Figure 4.3 plots the bootstrap MSEs, described in Section 3.3, for the three
areas with highest number of fires (Viana 1, Terra de Tribes and Viana 2). The
number of bootstrap replicates is $B=500$. The mean of the MSEs for the three
areas is $3.978$ in Model 0 and $1.219$ in Model ST1. Then, a clear
improvement is achieved when one uses the area-level spatio-temporal Poisson
mixed model ST1.

---
Figure 4.3: Bootstrap MSE estimates for the three areas with highest number of
fires.

##  5 Application to poverty data

Based on Model S1, as a particular case of Model ST1, this section applies the
developed statistical methodology to estimating women poverty proportions by
counties of Galicia during 2013. The data are taken from the Spanish living
conditions survey (SLCS), which provides annual information about the
household income received during the year prior to that of the interview.
Since there is only one time period (year 2013), this section removes the $t$
index.

For every individual, the equivalent personal income is obtained by dividing
the annual household net income by the equivalent total of household members,
which is obtained as a weighted sum. The poverty line is defined as a
percentage (currently Eurostat fixed it in 60%) of the median of the
equivalent personal incomes in the whole country. A person is defined as poor
(poverty=1) if his/her equivalent personal income is lower than the poverty
line and is defined as not poor (poverty=0) otherwise. The poverty proportion
of a territory is the average of the poverty variable. If the average is
restricted to women, then the corresponding poverty proportion for women is
obtained.

Spain is hierarchically divided in autonomous communities, provinces, counties
and municipalities. The SLCS is designed to obtain reliable statistics for
autonomous communities. Therefore, the usual direct estimators are not precise
enough for studying the poverty at a more disaggregated level, like provinces
or counties (small areas). We use the introduced SAE methodology for
constructing model-based predictors that give more accurate estimations of
women poverty proportions for each county area and find differences between
them.

In Galicia there are 53 counties, but in four of them there are no available
data. Therefore, the number of considered domains is $D=49$. The performance
of Model S1 depends on the choice of the proximity matrix $W$. Three different
choices are tested: common borders, distances and $k$-nearest neighbours. In
the first option (common border), two domains are neighbours if they have a
common delimitation. The last two options consider the Euclidean distance
between the centroids of the counties. The second option sets up a proximity
measure by taking the inverse of the distance between domains. The last option
applies $k$-nearest neighbours with $k=2$ and $3$. After analysing the
different possibilities in terms of model diagnostics, the first option is
selected because it is the one giving the best results. Figure 5.1 shows the
proximity map that determines the proximity matrix $\hbox{\boldmath$W$}^{0}$,
i.e. it provides for each domain, which are its neighbours. See Section 2 for
more details on the construction of the proximity matrix
$\hbox{\boldmath$W$}^{0}$ and $W$.

---
Figure 5.1: Proximity map for each domain $d$ ($d=1,\ldots,D$).

The target variable, $y_{d}$, counts the number of women under the poverty
line in the domain $d$ and $\nu_{d}$ is the corresponding sample size of
women. The minimum, median and maximum values of $\nu_{d}$ are $19$, $152$ and
$1384$, respectively. The minimum has been reached in the south east of the
region, while the median and maximum belong to the south west. The auxiliary
variables are the proportions of women who are unemployed ($lab2$) and who
have completed university education ($edu3$). These variables are calculated
by averaging direct estimates of the four Spanish employment surveys carried
out in 2013, all of them having bigger sample sizes than the SLCS. We first
fit Model 1 to the data $(y_{d},\nu_{d},\hbox{\boldmath$x$}_{d})$,
$d=1,\ldots,D$ and we apply Moran’s I test for spatial autocorrelation. As the
obtained $p$-value is lower than $0.001$, we assume that $y_{d}$,
$d=1,\ldots,D$, follows Model S1 and we fit this model to the data.

Table 5.1 provides a descriptive analysis for the logarithm of the response
variable, $\log y$, and the considered covariates. Specifically, it presents
their mean, standard deviation (sd), median, minimum and maximum values (min
and max) and the correlation (cor) between the covariates and $\log y$.

Table 5.1: Descriptive analysis. Variable | mean | sd | median | min | max | cor
---|---|---|---|---|---|---
$\log y$ | 3.13 | 1.15 | 3.18 | 0.00 | 5.51 | $-$
$lab2$ | 0.10 | 0.04 | 0.10 | 0.02 | 0.21 | 0.32
$edu3$ | 0.15 | 0.07 | 0.14 | 0.03 | 0.33 | 0.42

Table 5.2 presents the significant estimates ($p$-value $<0.05$) of the fixed
regression parameters under Model S1 and their standard errors, $z$-values and
$p$-values. The spatial autocorrelation parameter is estimated by applying
Moran’s I statistic (2.6) over the Pearson residuals of Model 0 and the
remaining model parameters are given as a solution of the system formed by the
first $p+1$ MM equations in (2.3).

Table 5.2: MM estimates of regression parameters under Model S1. Variable | Estimate | s.e. |  $z$-value |  $p$-value
---|---|---|---|---
$Intercept$ | -1.8803 | 0.1515 | -12.4086 |  $<$ 0.001
$lab2$ | 2.9848 | 1.2097 | 2.4689 | 0.0136
$edu3$ | -1.3809 | 0.5033 | -2.7445 | 0.0061

Taking into account the signs of the estimates, the auxiliary variable $lab2$
(proportion of unemployed women), is directly related to the response variable
while $edu3$ (proportion of women with university level of education), helps
to decrease the counts of women under the poverty line. Each domain $d$,
$d=1,\ldots,D$, has a random intercept with distribution $N(0,\phi^{2})$,
where $\hat{\phi}=0.130$. The $95\%$ percentile bootstrap confidence interval
for the standard deviation parameter is $(0.001,0.331)$. The estimated
autocorrelation parameter is $\hat{\rho}=0.324$. To test the null hypothesis
$H_{0}:\phi^{2}=0$, Algorithm 1 of Section 2 is applied. The obtained
$p$-value is $0.018$. Then, the null hypothesis is rejected at the level
$\alpha=0.05$. The Algorithm 2 of Section 2 adapted to case $T=1$, and which
compares Model S1 vs. Model 1, is applied to test $H_{0}:\rho=0$. The obtained
bootstrap $p$-value is $0.001$. Taking $\alpha=0.05$, the bootstrap test
concludes that the autocorrelation parameter $\rho$ is significantly different
from 0. Therefore, it recommends fitting Model S1 to the data, instead of
Model 1.

We compare the performance of Model S1 and Model 0, with only fixed-effects
(or equivalently, with $\phi=0$). For this sake, we fit Model 0 to the same
data as Model S1. Figure 5.2 plots the Pearson residuals

$r_{d}^{P}=\frac{y_{d}-\nu_{d}\hat{p}_{d}}{\sqrt{\nu_{d}\hat{p}_{d}}},\quad
d=1,\ldots,D,$

of the synthetic estimator
$\hat{p}_{d}^{syn}(\hat{\hbox{\boldmath$\beta$}})=\exp\\{\hbox{\boldmath$x$}_{d}\hat{\hbox{\boldmath$\beta$}}\\}$
based on Model 0 (left), and of the EBP approximation
$\hat{p}_{d}^{a}(\hat{\hbox{\boldmath$\theta$}})$ based on Model S1 (right).
In both cases, the distribution of the Pearson residuals is symmetrical around
0. In addition, the plots suggest a clear improvement when one uses an area-
level Poisson mixed model that incorporates SAR(1) domain effects, since it is
able to better capture the underlying spatial correlation structure. Figure
5.2 shows that most of the residuals of Model S1 take values in the interval
$(-2,2)$, while the residuals of Model 0 take values mainly in the interval
$(-3,3)$. In addition the tested hypotheses on $\rho$ and $\phi$, the
conclusion is again that Model S1 is more appropriated to fit the women
poverty data in Galicia by counties in $2013$.

|
---|---
Figure 5.2: Pearson residuals of the synthetic estimator based on Model 0
(left) and of the EBP approximation based on Model S1 (right).

Figure 5.3 (left) compares the behaviour of the EBP and direct estimations.
The direct estimators are calculated by using the Hájeck formula with the
officially calibrated sampling weights. The domains are sorted by the sample
sizes $\nu_{d}$’s. The direct estimators present oscillations of large
amplitude, while the EBPs have a smoother behaviour, which is something
preferred by the statistical offices when publishing estimations. As the
sample size increases, both sets of estimates tend to overlap.

Figure 5.3 (right) plots the relative root-MSEs (RRMSE) of the EBPs based on
Model S1 and the relative root-variances of the direct estimators. The RRMSEs
of the EBPs are estimated by using the bootstrap procedure of Section 3.3 with
$B=500$ replicates. For the direct estimators, we apply the formula (10.3.8)
of Särndal et al. (1992). The direct estimates have high variability,
specially for small sample sizes. As above, when the sample size $\nu_{d}$
increases, both accuracy measures follow the same pattern. The averages of the
relative root-variances of the direct estimator and of the RRMSEs of the EBP
are $0.2595$ and $0.1323$, respectively. According to these results, we
conclude that the EBP performs better.

|
---|---
Figure 5.3: Direct estimates and EBPs of poverty proportions $p_{d}$ (left)
and relative root-MSEs (right) for women in $2013$.

Figure 5.4 (left) maps the EBPs of $p_{d}$ for women in $2013$. The regions
where there is no data, are in white. Model S1 gives the following predictions
of women poverty proportions: 1 county with poverty proportion $p_{d}\leq
0.12$, 12 counties with $0.12<p_{d}\leq 0.15$, 24 counties with
$0.15<p_{d}\leq 0.18$ and 12 counties with $p_{d}>0.18$. Highest levels of
poverty are found in the south and west of the community. On the other hand,
the counties with the lowest estimated poverty proportion are located in the
north-east of the region.

|
---|---
Figure 5.4: Poverty proportion EBPs for women based on Model S1 (left) and
RRMSEs (right) in Galicia during $2013$.

Figure 5.4 (right) maps the RRMSE estimates of the EBP of $p_{d}$ by counties
in $2013$. There are 8 counties with RRMSE $\leq 10\%$, 8 counties with
$10\%<$ RRMSE $\leq 13\%$, 19 counties with $13\%<$ RRMSE $\leq 16\%$ and 14
counties with RRMSE $>16\%$. The highest values are found in the north-east of
the region. Their minimum and maximum are $8.82\%$ and $18.49\%$,
respectively. As the highest RRMSE is lower than $20\%$, these estimates could
be accepted for publication by statistical offices.

##  6 Concluding remarks

This paper introduces Model ST1, which is an area-level spatio-temporal
Poisson mixed model with SAR(1)-correlated domain random effects and
independent domain-time random effects. It contains the area-level Poisson
mixed model with SAR(1) domain effects (Model S1) and the area-level Poisson
mixed model with independent time effects (Model T1) as particular or limit
cases. The MM algorithm is obtained to fit the model parameters. It is based
on the method of simulated moments proposed proposed by Jiang (1998).

Two predictors of the target parameter $p_{dt}$ are proposed. They are the EBP
and the plug-in predictor. The EBP is asymptotically unbiased with minimum MSE
and it is computationally faster than the plug-in predictor, since the last
one requires calculating two EBPs ($\hat{\hbox{\boldmath$v$}}_{1}$ and
$\hat{\hbox{\boldmath$v$}}_{2}$). Further, an approximation to the EBP is
proposed to avoid the computational burden.

The developed methodology is applied to the data set of wildfires in Galicia
by forest areas and months during 2007-2008. Because of their better
properties, the EBP is employed. This predictor is also compared against the
synthetic estimator obtained under Model 0, with only fixed effects. A clear
improvement is achieved when one uses a more complex model incorporating
random effects. In addition, this paper recommends using Model ST1 to analyse
wildfires in Galicia since the bootstrap test of Algorithm 2 yielded
$p$-values close to 0.

The conclusion is that forest fires tend to be concentrated in coastal areas
and in the south of the region. An important increase in wildfires is observed
in September and October $2007$. The introduced bootstrap MSE estimator is
considered as accuracy measure of the proposed EBP. It is achieved a clear
improvement when using the proposed model and estimator against the classical
Poisson regression model.

We also use the EBP approximation for estimating women poverty proportions in
Galician counties. The data are taken from the 2013 SLCS. As both spatial
correlation tests, the Moran’s I test and the modified Algorithm 2, indicates
spatial correlation, we fit Model S1 to the data. In addition, the proposed
model-based predictor is compared against the direct estimator. The EBPs of
the women poverty proportion are smoother. As the RRMSE of the direct
estimator is too high when the sample size $\nu_{d}$ is small, it is
preferable to use the EBP. The estimates based on Model S1 suggests that the
highest levels of women poverty are found in the south and west of the region.
The average percentage of women poverty is $16.89\%$ and its average error is
$13.23\%$.

It is worth mentioning the programming problems of the predictors constructed
under the introduced spatio-temporal models, which leads us to use numerical
approaches that introduce an additional source of error. However, if there is
a spatial correlation, the approximate EBPs behave better than other
predictors that do not take into account this information.

##  Appendix A Appendix: MM equations

First, we recall that the moment generation function of $Y\sim
N(\mu,\sigma^{2})$ is

$\Psi(t;\mu,\sigma^{2})=\mathbb{E}\left[e^{tY}\right]=\exp\big{\\{}\mu
t+\frac{1}{2}\,\sigma^{2}t^{2}\big{\\}}.$

For ease of exposition, the elements of $\Gamma$ and its derivatives are
denoted by $\gamma_{d_{1}d_{2}}=\gamma_{d_{1}d_{2}}(\rho)$ and
$\dot{\gamma}_{d_{1}d_{2}}=\dot{\gamma}_{d_{1}d_{2}}(\rho)$, respectively. The
calculations start with the first $p$ MM equations. The expectation of
$y_{dt}$ is

$\displaystyle\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}]$
$\displaystyle=\mathbb{E}_{v}\big{[}E_{\hbox{\boldmath$\theta$}}[y_{dt}|\hbox{\boldmath$v$}]\big{]}=\mathbb{E}_{v}[\nu_{dt}p_{dt}]=\mathbb{E}_{v}\left[\nu_{dt}\exp\left\\{\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\right\\}\right]$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\nu_{dt}\exp\left\\{\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\right\\}\,f(v_{1,d})f(v_{2,dt})\,dv_{1,d}dv_{2,dt}$
$\displaystyle=\int_{-\infty}^{\infty}\nu_{dt}\exp\Big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}\phi_{2}^{2}+\phi_{1}v_{1,d}\Big{\\}}f_{v}(v_{1,d})\,dv_{1,d}$
$\displaystyle=\nu_{dt}\exp\Big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}\big{(}\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2}\big{)}\Big{\\}}.$

Therefore, the first $p$ MM equations are

$\displaystyle f_{k}(\hbox{\boldmath$\theta$})$
$\displaystyle=\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}\nu_{dt}\exp\Big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}\big{(}\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2}\big{)}\Big{\\}}x_{dtk}$
$\displaystyle-\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}y_{dt}x_{dtk},\,\,k=1,\ldots,p.$

The derivatives of $\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}]$ are

$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}]}{\partial\beta_{k}}$
$\displaystyle=\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}x_{dtk},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}]}{\partial\phi_{1}}$
$\displaystyle=\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{1}\gamma_{dd},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}]}{\partial\phi_{2}}$
$\displaystyle=\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{2},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}]}{\partial\rho}$
$\displaystyle=\frac{1}{2}\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{1}^{2}\dot{\gamma}_{dd}.$

The expectation of $y_{dt}^{2}$ is
$\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}]=\mathbb{E}_{v}\big{[}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}|\hbox{\boldmath$v$}]\big{]}$,
where

$\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}|\hbox{\boldmath$v$}]=\mbox{var}_{\hbox{\boldmath$\theta$}}[y_{dt}|\hbox{\boldmath$v$}]+\mathbb{E}_{\hbox{\boldmath$\theta$}}^{2}[y_{dt}|\hbox{\boldmath$v$}]=\nu_{dt}p_{dt}+\nu_{dt}^{2}p_{dt}^{2}.$

Therefore

$\displaystyle\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}]=\mathbb{E}_{v}\big{[}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}|\hbox{\boldmath$v$}]\big{]}$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\nu_{dt}p_{dt}f(v_{1,d})f(v_{2,dt})\,dv_{1,d}dv_{2,dt}$
$\displaystyle+\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\nu_{dt}^{2}p_{dt}^{2}f(v_{1,d})f(v_{2,dt})\,dv_{1,d}dv_{2,dt}=S_{1}+S_{2},$

where

$\displaystyle S_{2}$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\nu_{dt}^{2}p_{dt}^{2}f(v_{2,dt})f(v_{1,d})\,dv_{2,dt}dv_{1,d}$
$\displaystyle=\nu_{dt}^{2}\int_{-\infty}^{\infty}\left[\int_{-\infty}^{\infty}\exp\big{\\{}2\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+2\phi_{1}v_{1,d}+2\phi_{2}v_{2,dt}\big{\\}}f(v_{2,dt})\,dv_{2,dt}\right]f(v_{1,d})\,dv_{1,d}$
$\displaystyle=\nu_{dt}^{2}\int_{-\infty}^{\infty}\exp\big{\\{}2(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{2}^{2})+2\phi_{1}v_{1,d}\big{\\}}f(v_{1,d})\,dv_{1,d}$
$\displaystyle=\nu_{dt}^{2}\exp\big{\\{}2(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}.$

Then, the expectation $\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}]$ is

$\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}]=\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}+\nu_{dt}^{2}\exp\big{\\{}2(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}},$

and as a consequence, the $(p+2)$-th MM equation is

$\displaystyle f_{p+2}(\hbox{\boldmath$\theta$})$
$\displaystyle=\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}\bigg{\\{}\nu_{dt}\exp\Big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\Big{\\}}+\nu_{dt}^{2}\exp\big{\\{}2(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\bigg{\\}}$
$\displaystyle-\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}y_{dt}^{2}.$

The derivatives of $\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}]$ are

$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}]}{\partial\beta_{k}}$
$\displaystyle=\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}x_{dtk}+2\nu_{dt}^{2}\exp\big{\\{}2(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}x_{dtk},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d}^{2}]}{\partial\phi_{1}}$
$\displaystyle=\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{1}\gamma_{dd}+4\nu_{dt}^{2}\exp\big{\\{}2(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{1}\gamma_{dd},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d}^{2}]}{\partial\phi_{2}}$
$\displaystyle=\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{2}+4\nu_{dt}^{2}\exp\big{\\{}2(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{2},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d}^{2}]}{\partial\rho}$
$\displaystyle=\frac{1}{2}\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{1}^{2}\dot{\gamma}_{dd}+2\nu_{dt}^{2}\exp\big{\\{}2(\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\phi_{1}^{2}\dot{\gamma}_{dd}.$

The expectation of $y_{d.}^{2}$ is
$\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]=\mathbb{E}_{v}\big{[}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}|\hbox{\boldmath$v$}]\big{]}$,
where

$y_{d.}^{2}=\sum_{t=1}^{T}y_{dt}^{2}+\sum_{t_{1}\neq
t_{2}}y_{dt_{1}}y_{dt_{2}},\quad\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}|\hbox{\boldmath$v$}]=\mbox{var}_{\hbox{\boldmath$\theta$}}[y_{dt}|\hbox{\boldmath$v$}]+\mathbb{E}_{\hbox{\boldmath$\theta$}}^{2}[y_{dt}|\hbox{\boldmath$v$}]=\nu_{dt}p_{dt}+\nu_{dt}^{2}p_{dt}^{2}.$

The expectation of $y_{d.}^{2}$, conditionally on the random effects $v$, is

$\displaystyle\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}|\hbox{\boldmath$v$}]$
$\displaystyle=\sum_{t=1}^{T}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}|\hbox{\boldmath$v$}]+\sum_{t_{1}\neq
t_{2}}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt_{1}}|\hbox{\boldmath$v$}]\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt_{2}}|\hbox{\boldmath$v$}]$
$\displaystyle=\sum_{t=1}^{T}\big{\\{}\nu_{dt}p_{dt}+\nu_{dt}^{2}p_{dt}^{2}\big{\\}}+\sum_{t_{1}\neq
t_{2}}\nu_{dt_{1}}p_{dt_{1}}\nu_{dt_{2}}p_{dt_{2}}.$

Therefore

$\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]=\sum_{t=1}^{T}\nu_{dt}\mathbb{E}_{v}[p_{dt}]+\sum_{t=1}^{T}\nu_{dt}^{2}\mathbb{E}_{v}[p_{dt}^{2}]+\sum_{t_{1}\neq
t_{2}}\nu_{dt_{1}}\nu_{dt_{2}}\mathbb{E}_{v}[p_{dt_{1}}p_{dt_{2}}],$

where the expectation of $p_{dt_{1}}p_{dt_{2}}$ is

$\displaystyle\mathbb{E}_{v}[p_{dt_{1}}p_{dt_{2}}]$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\exp\\{(\hbox{\boldmath$x$}_{dt_{1}}+\hbox{\boldmath$x$}_{dt_{2}})\hbox{\boldmath$\beta$}+2\phi_{1}v_{1,d}+\phi_{2}v_{2,dt_{1}}+\phi_{2}v_{2,dt_{2}}\\}$
$\displaystyle\cdot
f(v_{2,dt_{1}})f(v_{2,dt_{2}})f(v_{1,d})\,dv_{2,dt_{1}}dv_{2,dt_{2}}dv_{1,d}$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\exp\\{(\hbox{\boldmath$x$}_{dt_{1}}+\hbox{\boldmath$x$}_{dt_{2}})\hbox{\boldmath$\beta$}+\frac{1}{2}\phi_{2}^{2}+\phi_{2}v_{2,dt_{1}}+2\phi_{1}v_{1,d}\\}$
$\displaystyle\cdot f(v_{2,dt_{1}})f(v_{1,d})\,dv_{2,dt_{1}}dv_{1,d}$
$\displaystyle=\int_{-\infty}^{\infty}\exp\\{(\hbox{\boldmath$x$}_{dt_{1}}+\hbox{\boldmath$x$}_{dt_{2}})\hbox{\boldmath$\beta$}+\frac{1}{2}\phi_{2}^{2}+\frac{1}{2}\phi_{2}^{2}+2\phi_{1}v_{1,d}\\}f(v_{1,d})\,dv_{1,d}$
$\displaystyle=\exp\\{(\hbox{\boldmath$x$}_{dt_{1}}+\hbox{\boldmath$x$}_{dt_{2}})\hbox{\boldmath$\beta$}+\phi_{2}^{2}+2\phi_{1}^{2}\gamma_{dd}\\}.$

Then, the expectation of $y_{d.}^{2}$ is

$\displaystyle\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]$
$\displaystyle=\sum_{t=1}^{T}\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}+\sum_{t=1}^{T}\nu_{dt}^{2}\exp\big{\\{}2\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+2(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}$
$\displaystyle+\sum_{t_{1}\neq
t_{2}}\nu_{dt_{1}}\nu_{dt_{2}}\exp\\{(\hbox{\boldmath$x$}_{dt_{1}}+\hbox{\boldmath$x$}_{dt_{2}})\hbox{\boldmath$\beta$}+2\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2}\\}$
$\displaystyle\pm\sum_{t=1}^{T}\nu_{dt}^{2}\exp\big{\\{}2\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+2\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2}\big{\\}}$
$\displaystyle=\sum_{t=1}^{T}\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}+\sum_{t=1}^{T}\nu_{dt}^{2}\exp\big{\\{}2\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+2(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}$
$\displaystyle-\sum_{t=1}^{T}\nu_{dt}^{2}\exp\big{\\{}2\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+2\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2}\big{\\}}+\left(\sum_{t=1}^{T}\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\frac{1}{2}\phi_{2}^{2}\big{\\}}\right)^{2},$

and as a consequence, the $(p+1)$-th MM equation is

$\displaystyle f_{p+1}(\hbox{\boldmath$\theta$})$
$\displaystyle=\frac{1}{D}\sum_{d=1}^{D}\left\\{\vphantom{\left(\sum_{t=1}^{T}\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\frac{1}{2}\phi_{2}^{2}\big{\\}}\right)^{2}}\sum_{t=1}^{T}\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\right.$
$\displaystyle+\left.\big{(}e^{\phi_{2}^{2}}-1\big{)}\sum_{t=1}^{T}\nu_{dt}^{2}\exp\big{\\{}2\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+2\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}}\right.$
$\displaystyle+\left.\left(\sum_{t=1}^{T}\nu_{dt}\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\frac{1}{2}\phi_{2}^{2}\big{\\}}\right)^{2}\right\\}-\frac{1}{D}\sum_{d=1}^{D}y_{d.}^{2}.$

The derivatives of $\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]$ are

$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]}{\partial\beta_{k}}$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\nu_{dt}P_{dt}x_{dtk}+2\sum_{t=1}^{T}\nu_{dt}^{2}Q_{dt}x_{dtk}-2\sum_{t=1}^{T}\nu_{dt}^{2}R_{dt}x_{dtk}$
$\displaystyle+$ $\displaystyle
2\bigg{(}\sum_{t=1}^{T}\nu_{dt}S_{dt}\bigg{)}\sum_{t=1}^{T}\nu_{dt}S_{dt}x_{dtk},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]}{\partial\phi_{1}}$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\nu_{dt}P_{dt}\phi_{1}\gamma_{dd}+4\sum_{t=1}^{T}\nu_{dt}^{2}Q_{dt}\phi_{1}\gamma_{dd}-4\sum_{t=1}^{T}\nu_{dt}^{2}R_{dt}\phi_{1}\gamma_{dd}$
$\displaystyle+$ $\displaystyle
4\bigg{(}\sum_{t=1}^{T}\nu_{dt}S_{dt}\bigg{)}\sum_{t=1}^{T}\nu_{dt}S_{dt}\phi_{1}\gamma_{dd},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]}{\partial\phi_{2}}$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\nu_{dt}P_{dt}\phi_{2}+4\sum_{t=1}^{T}\nu_{dt}^{2}Q_{dt}\phi_{2}-2\sum_{t=1}^{T}\nu_{dt}^{2}R_{dt}\phi_{2}$
$\displaystyle+$ $\displaystyle
2\bigg{(}\sum_{t=1}^{T}\nu_{dt}S_{dt}\bigg{)}\sum_{t=1}^{T}\nu_{dt}S_{dt}\phi_{2},$
$\displaystyle\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]}{\partial\rho}$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\frac{1}{2}\nu_{dt}P_{dt}\phi_{1}^{2}\dot{\gamma}_{dd}+2\sum_{t=1}^{T}\nu_{dt}^{2}Q_{dt}\phi_{1}^{2}\dot{\gamma}_{dd}-2\sum_{t=1}^{T}\nu_{dt}^{2}R_{dt}\phi_{1}^{2}\dot{\gamma}_{dd}$
$\displaystyle+$ $\displaystyle
2\bigg{(}\sum_{t=1}^{T}\nu_{dt}S_{dt}\bigg{)}\sum_{t=1}^{T}\nu_{dt}S_{dt}\phi_{1}^{2}\dot{\gamma}_{dd},$

where

$\displaystyle P_{dt}$ $\displaystyle=$
$\displaystyle\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\frac{1}{2}(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}},\quad
Q_{dt}=\exp\big{\\{}2\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+2(\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2})\big{\\}},$
$\displaystyle R_{dt}$ $\displaystyle=$
$\displaystyle\exp\big{\\{}2\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+2\phi_{1}^{2}\gamma_{dd}+\phi_{2}^{2}\big{\\}},\quad
S_{dt}=\exp\big{\\{}\hbox{\boldmath$x$}_{dt}\hbox{\boldmath$\beta$}+\phi_{1}^{2}\gamma_{dd}+\frac{1}{2}\phi_{2}^{2}\big{\\}}.$

The expectation of $y_{d_{1}.}y_{d_{2}.}$ is

$\displaystyle\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d_{1}.}y_{d_{2}.}]$
$\displaystyle=\mathbb{E}_{v}\big{[}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d_{1}.}y_{d_{2}.}|\hbox{\boldmath$v$}]\big{]}=\mathbb{E}_{v}\big{[}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d_{1}.}|v_{1,d_{1}},\hbox{\boldmath$v$}_{2,d_{1}}]\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d_{2}.}|v_{1,d_{2}},\hbox{\boldmath$v$}_{2,d_{2}}]\big{]}$
$\displaystyle=\sum_{t_{1}=1}^{T}\sum_{t_{2}=1}^{T}\mathbb{E}_{v}\big{[}\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d_{1}t_{1}}|v_{1,d_{1}},v_{2,d_{1}t_{1}}]\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d_{2}t_{2}}|v_{1,d_{2}},v_{2,d_{2}t_{2}}]\big{]}$
$\displaystyle=\sum_{t_{1}=1}^{T}\sum_{t_{2}=1}^{T}\nu_{d_{1}t_{1}}\nu_{d_{2}t_{2}}\mathbb{E}_{v}\big{[}p_{d_{1}t_{1}}p_{d_{2}t_{2}}\big{]}.$

By defining
$\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})=\mathbb{E}_{v}\big{[}p_{d_{1}t_{1}}p_{d_{2}t_{2}}\big{]}$,
it holds that

$\displaystyle\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{4}}\exp\big{\\{}(\hbox{\boldmath$x$}_{d_{1}t_{1}}+\hbox{\boldmath$x$}_{d_{2}t_{2}})\hbox{\boldmath$\beta$}+\phi_{1}(v_{1,d_{1}}+v_{1,d_{2}})+\phi_{2}(v_{2,d_{1}t_{1}}+v_{2,d_{2}t_{2}})\big{\\}}$
$\displaystyle\cdot
f(v_{2,d_{2}t_{2}})dv_{2,d_{2}t_{2}}\,f(v_{2,d_{1}t_{1}})dv_{2,d_{1}t_{1}}f(v_{1,d_{2}}|v_{1,d_{1}})dv_{1,d_{2}}\,f(v_{1,d_{1}})\,dv_{1,d_{1}}$
$\displaystyle=\int_{\mathbb{R}^{3}}\exp\big{\\{}(\hbox{\boldmath$x$}_{d_{1}t_{1}}+\hbox{\boldmath$x$}_{d_{2}t_{2}})\hbox{\boldmath$\beta$}+\phi_{1}(v_{1,d_{1}}+v_{1,d_{2}})+\phi_{2}v_{2,d_{1}t_{1}}+\frac{1}{2}\phi_{2}^{2}\big{\\}}$
$\displaystyle\cdot
f(v_{2,d_{1}t_{1}})dv_{2,d_{1}t_{1}}\,f(v_{1,d_{2}}|v_{1,d_{1}})dv_{1,d_{2}}\,f(v_{1,d_{1}})\,dv_{1,d_{1}}.$

Therefore,

$\displaystyle\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})$
$\displaystyle=\int_{\mathbb{R}^{2}}\exp\big{\\{}(\hbox{\boldmath$x$}_{d_{1}t_{1}}+\hbox{\boldmath$x$}_{d_{2}t_{2}})\hbox{\boldmath$\beta$}+\phi_{1}(v_{1,d_{1}}+v_{1,d_{2}})+\frac{1}{2}\phi_{2}^{2}+\frac{1}{2}\phi_{2}^{2}\big{\\}}$
$\displaystyle\cdot
f(v_{1,d_{2}}|v_{1,d_{1}})dv_{1,d_{2}}\,f(v_{1,d_{1}})\,dv_{1,d_{1}}$
$\displaystyle=\int_{\mathbb{R}}\exp\Big{\\{}(\hbox{\boldmath$x$}_{d_{1}t_{1}}+\hbox{\boldmath$x$}_{d_{2}t_{2}})\hbox{\boldmath$\beta$}+\phi_{1}v_{1,d_{1}}+\frac{\gamma_{d_{1}d_{2}}}{\gamma_{d_{1}d_{1}}}v_{1,d_{1}}\phi_{1}$
$\displaystyle+\frac{1}{2}\Big{(}\gamma_{d_{2}d_{2}}-\frac{\gamma_{d_{1}d_{2}}^{2}}{\gamma_{d_{1}d_{1}}}\Big{)}\phi_{1}^{2}+\phi_{2}^{2}\Big{\\}}f(v_{1,d_{1}})\,dv_{1,d_{1}}$
$\displaystyle=\exp\Big{\\{}(\hbox{\boldmath$x$}_{d_{1}t_{1}}+\hbox{\boldmath$x$}_{d_{2}t_{2}})\hbox{\boldmath$\beta$}+\frac{1}{2}\Big{(}1+\frac{\gamma_{d_{1}d_{2}}}{\gamma_{d_{1}d_{1}}}\Big{)}^{2}\gamma_{d_{1}d_{1}}\phi_{1}^{2}+\frac{1}{2}\Big{(}\gamma_{d_{2}d_{2}}-\frac{\gamma_{d_{1}d_{2}}^{2}}{\gamma_{d_{1}d_{1}}}\Big{)}\phi_{1}^{2}+\phi_{2}^{2}\Big{\\}}$
$\displaystyle=\exp\Big{\\{}(\hbox{\boldmath$x$}_{d_{1}t_{1}}+\hbox{\boldmath$x$}_{d_{2}t_{2}})\hbox{\boldmath$\beta$}+\frac{1}{2}\,\phi_{1}^{2}(\gamma_{d_{1}d_{1}}+2\gamma_{d_{1}d_{2}}+\gamma_{d_{2}d_{2}})+\phi_{2}^{2}\Big{\\}}.$

Therefore, the $(p+3)$-th MM equation is

$\displaystyle f_{p+3}(\hbox{\boldmath$\theta$})$
$\displaystyle=\frac{1}{D(D-1)}\sum_{d_{1}\neq
d_{2}}^{D}\sum_{t_{1}=1}^{T}\sum_{t_{2}=1}^{T}\nu_{d_{1}t_{1}}\nu_{d_{2}t_{2}}\exp\Big{\\{}(\hbox{\boldmath$x$}_{d_{1}t_{1}}+\hbox{\boldmath$x$}_{d_{2}t_{2}})\hbox{\boldmath$\beta$}$
$\displaystyle+\frac{1}{2}\,\phi_{1}^{2}(\gamma_{d_{1}d_{1}}+2\gamma_{d_{1}d_{2}}+\gamma_{d_{2}d_{2}})+\phi_{2}^{2}\Big{\\}}-\frac{1}{D(D-1)}\sum_{d_{1}\neq
d_{2}}^{D}\sum_{t_{1}=1}^{T}\sum_{t_{2}=1}^{T}\nu_{d_{1}t_{1}}y_{d_{1}.}y_{d_{2}.}$

The derivatives of
$\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})$ are

$\displaystyle\frac{\partial\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})}{\partial\beta_{k}}$
$\displaystyle=\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})(x_{d_{1}t_{1}k}+x_{d_{2}t_{2}k}),$
$\displaystyle\frac{\partial\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})}{\partial\phi_{1}}$
$\displaystyle=\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})\phi_{1}(\gamma_{d_{1}d_{1}}+\gamma_{d_{2}d_{2}}+2\gamma_{d_{1}d_{2}}),$
$\displaystyle\frac{\partial\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})}{\partial\phi_{2}}$
$\displaystyle=2\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})\phi_{2},$
$\displaystyle\frac{\partial\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})}{\partial\rho}$
$\displaystyle=\frac{1}{2}\varphi_{d_{1},d_{2}}(\hbox{\boldmath$\theta$})\phi_{1}^{2}(\dot{\gamma}_{d_{1}d_{1}}+\dot{\gamma}_{d_{2}d_{2}}+2\dot{\gamma}_{d_{1}d_{2}}).$

The elements of the Jacobian matrix are

$\displaystyle H_{kr}$ $\displaystyle=\frac{\partial
f_{k}(\hbox{\boldmath$\theta$})}{\partial\theta_{r}}=\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}]}{\partial\theta_{r}}\,x_{dtk},\quad
k=1,\ldots,p,\,r=1,\ldots,p+3,$ $\displaystyle H_{p+1r}$
$\displaystyle=\frac{\partial
f_{p+1}(\hbox{\boldmath$\theta$})}{\partial\theta_{r}}=\frac{1}{D}\sum_{d=1}^{D}\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{d.}^{2}]}{\partial\theta_{r}},\quad
r=1,\ldots,p+3,$ $\displaystyle H_{p+2r}$ $\displaystyle=\frac{\partial
f_{p+2}(\hbox{\boldmath$\theta$})}{\partial\theta_{r}}=\frac{1}{DT}\sum_{d=1}^{D}\sum_{t=1}^{T}\frac{\partial\mathbb{E}_{\hbox{\boldmath$\theta$}}[y_{dt}^{2}]}{\partial\theta_{r}},\quad
r=1,\ldots,p+3,$ $\displaystyle H_{p+3r}$ $\displaystyle=\frac{\partial
f_{p+3}(\hbox{\boldmath$\theta$})}{\partial\theta_{r}}=\frac{1}{D(D-1)}\sum_{d_{1}\neq
d_{2}}^{D}\sum_{t_{1}=1}^{T}\sum_{t_{2}=1}^{T}\nu_{d_{1}t_{1}}\nu_{d_{2}t_{2}}\frac{\partial\varphi_{d_{1}d_{2}}^{t_{1}t_{2}}(\hbox{\boldmath$\theta$})}{\partial\theta_{r}},\quad
r=1,\ldots,p+3.$

##  Appendix B Appendix: Simulations

This section presents two model-based simulation experiments. The first one
studies the the behaviour of the MM fitting algorithm while the second one
compares the performance of the two introduced predictors, i.e. the EBP and
the plug-in. The response variables are generated independently as
$y_{dt}|v_{1,d},v_{2,dt}\sim\mbox{Poisson}(\nu_{dt}p_{dt})$, where

$p_{dt}=\exp\\{\beta_{0}+x_{dt}\beta_{1}+\phi_{1}v_{1,d}+\phi_{2}v_{2,dt}\\},\,x_{dt}=\frac{d+t/T}{D},\,d=1,\ldots,D,\,t=1,\ldots,T.$

The domain random effects, $v_{1,d}\,(d=1,\ldots,D)$, are generated according
to a SAR(1) process, i.e.

$\hbox{\boldmath$v$}_{1}=\underset{1\leq d\leq
D}{\mbox{col}}(v_{1,d})=(\hbox{\boldmath$I$}_{D}-\rho\hbox{\boldmath$W$})^{-1}\hbox{\boldmath$u$}_{1},$

where $\hbox{\boldmath$I$}_{D}$ denotes the $D\times D$ identity matrix,
$\rho$ is the autocorrelation parameter,
$\hbox{\boldmath$W$}=(\omega_{ij})_{i,j=1,\ldots,D}$ is a proximity matrix and
$\hbox{\boldmath$u$}_{1}\sim N(\hbox{\boldmath$0$},\hbox{\boldmath$I$}_{D})$.
For the $D\times D$ proximity matrix $W$, a $7$-diagonal matrix is considered.
Let $k$ be the number of diagonals of $W$, then the number of upper and lower
diagonals is $m=\lfloor k/2\rfloor$, where $\lfloor k/2\rfloor$ denotes
integer part of $k/2$. The diagonals are denoted by $1$ (main) and $j$ (upper
and lower), $j=2,\ldots,m+1$. The diagonals are constructed in the following
way.

  * •

Diagonal 1 (main): if $|i-j|=0$, then $\omega_{ij}=0$.

  * •

Diagonal 2 (upper and lower): if $|i-j|=1$, then
$\omega_{ij}=\frac{1}{2}-\big{(}\frac{1}{2^{3}}+\frac{1}{2^{4}}\big{)}$.

  * •

Diagonals $3-(m+1)$ (upper and lower): if $|i-j|\in\\{2,\ldots,m\\}$, then
$\omega_{ij}=1/2^{|i-j|+1}$.

Then, the elements of diagonals 1, 2, 3 and 4 (upper and lower) are $0$,
$5/16$, $2/16$ and $1/16$ respectively. This rule does not apply to the first
and the last $m$ rows ($m=3$). For those rows, the numerators of the diagonal
elements are kept fixed ($0$, $5$, $2$ and $1$) and the denominators are
recalculated so that the sum of the row is $1$. Using this criterion, the
$9\times 9$ $7$-diagonal matrix $W$ is

$\hbox{\boldmath$W$}=\left(\begin{array}[]{ccccccccc}0&5/8&2/8&1/8&0&0&0&0&0\\\
5/13&0&5/13&2/13&1/13&0&0&0&0\\\ 2/15&5/15&0&5/15&2/15&1/15&0&0&0\\\ \hline\cr
1/16&2/16&5/16&0&5/16&2/16&1/16&0&0\\\ 0&1/16&2/16&5/16&0&5/16&2/16&1/16&0\\\
0&0&1/16&2/16&5/16&0&5/16&2/16&1/16\\\ \hline\cr
0&0&0&1/15&2/15&5/15&0&5/15&2/15\\\ 0&0&0&0&1/13&2/13&5/13&0&5/13\\\
0&0&0&0&0&1/8&2/8&5/8&0\\\ \end{array}\right).$

That $9\times 9$ $7$-diagonal matrix $W$ can be generalized to a $D\times D$
matrix by repeating $D-6$ times the weights of the central rows (i.e., $0$,
$5/16$, $2/16$ and $1/16$).

In both simulation experiments, we take $\beta_{0}=-3$, $\beta_{1}=0.8$,
$\phi_{1}=0.5$, $\phi_{2}=0.5$ and $\nu_{dt}=100$, $d=1,\ldots,D$,
$t=1,\ldots,T$. The simulation considers the scenarios $D=100$ and $T=4,8$ for
studying the influence of the time periods. For each scenario, it takes
$\rho=0.1,0.3,0.5$. This section runs the Monte Carlo simulation experiments
with $K=1000$ iterations.

###  B.1 Simulation 1

The target of Simulation 1 is to check the behaviour of the MM fitting
algorithm introduced in Section 2. Table B.1 presents the bias and Table B.2
the root mean squared error (RMSE) for the model parameters
$\theta\in\hbox{\boldmath$\theta$}=\\{\beta_{0},\beta_{1},\phi_{1},\phi_{2},\rho\\}$.
This section considers two options for estimating the vector of all model
parameters $\theta$. In the first option (Opt. 1), the vector
$\hat{\hbox{\boldmath$\theta$}}$ is obtained as a solution of the system of
$p+3$ nonlinear equations (2.3). In the second option (Opt. 2), $\hat{\rho}$
is given by calculating the Moran’s I measure over the predicted domain random
effects under Model T1.

Table B.1: Bias of the MM fitting algorithm under Model ST1. |  | $T=4$ |  | $T=8$
---|---|---|---|---
$\rho$ |  | Opt. 1 | Opt. 2 |  | Opt. 1 | Opt. 2
0.1 | $\hat{\beta}_{0}$ | 0.0115 | 0.0115 |  | 0.0171 | 0.0171
| $\hat{\beta}_{1}$ | -0.0145 | -0.0145 |  | -0.0233 | -0.0234
| $\hat{\phi}_{1}$ | -0.0219 | -0.0230 |  | -0.0240 | -0.0250
| $\hat{\phi}_{2}$ | -0.0098 | -0.0098 |  | -0.0066 | -0.0066
| $\hat{\rho}$ | -0.1689 | -0.0848 |  | -0.1593 | -0.0820
0.3 | $\hat{\beta}_{0}$ | 0.0196 | 0.0197 |  | 0.0150 | 0.0151
| $\hat{\beta}_{1}$ | -0.0225 | -0.0226 |  | -0.0168 | -0.0169
| $\hat{\phi}_{1}$ | -0.0144 | -0.0186 |  | -0.0129 | -0.0181
| $\hat{\phi}_{2}$ | -0.0099 | -0.0099 |  | -0.0087 | -0.0087
| $\hat{\rho}$ | -0.3236 | -0.1986 |  | -0.3441 | -0.1865
0.5 | $\hat{\beta}_{0}$ | 0.0285 | 0.0285 |  | 0.0394 | 0.0394
| $\hat{\beta}_{1}$ | -0.0387 | -0.0387 |  | -0.0666 | -0.0664
| $\hat{\phi}_{1}$ | 0.0101 | 0.0059 |  | 0.0062 | 0.0033
| $\hat{\phi}_{2}$ | -0.0123 | -0.0123 |  | -0.0111 | -0.0111
| $\hat{\rho}$ | -0.7090 | -0.2856 |  | -0.8400 | -0.2601

Both options behave similarly for the fixed effects and the variance
parameters. For these parameters, the variance is the most important term of
the MSE since bias is much smaller than RMSE. On the other hand, Opt. 2
produces more competitive estimates for the autocorrelation parameter $\rho$,
since it drastically reduces both bias and RMSE. For $\rho$, bias is the main
part of the MSE since it takes similar absolute values to the RMSE. Then, a
bias correction by bootstrap might be useful.

Table B.2: Root mean squared error of the MM fitting algorithm under Model ST1. |  | $T=4$ |  | $T=8$
---|---|---|---|---
$\rho$ |  | Opt. 1 | Opt. 2 |  | Opt. 1 | Opt. 2
0.1 | $\hat{\beta}_{0}$ | 0.1388 | 0.1387 |  | 0.1291 | 0.1291
| $\hat{\beta}_{1}$ | 0.2337 | 0.2337 |  | 0.2162 | 0.2161
| $\hat{\phi}_{1}$ | 0.0658 | 0.0659 |  | 0.0561 | 0.0563
| $\hat{\phi}_{2}$ | 0.0449 | 0.0449 |  | 0.0289 | 0.0289
| $\hat{\rho}$ | 0.1904 | 0.1129 |  | 0.1796 | 0.1090
0.3 | $\hat{\beta}_{0}$ | 0.1682 | 0.1682 |  | 0.1554 | 0.1554
| $\hat{\beta}_{1}$ | 0.2824 | 0.2823 |  | 0.2609 | 0.2609
| $\hat{\phi}_{1}$ | 0.0649 | 0.0647 |  | 0.0569 | 0.0545
| $\hat{\phi}_{2}$ | 0.0463 | 0.0463 |  | 0.0309 | 0.0309
| $\hat{\rho}$ | 0.3537 | 0.2144 |  | 0.3872 | 0.2039
0.5 | $\hat{\beta}_{0}$ | 0.2189 | 0.2189 |  | 0.2100 | 0.2100
| $\hat{\beta}_{1}$ | 0.3736 | 0.3736 |  | 0.3509 | 0.3509
| $\hat{\phi}_{1}$ | 0.0742 | 0.0650 |  | 0.0696 | 0.0530
| $\hat{\phi}_{2}$ | 0.0451 | 0.0451 |  | 0.0307 | 0.0307
| $\hat{\rho}$ | 0.7766 | 0.2998 |  | 0.9038 | 0.2755

###  B.2 Simulation 2

The second simulation experiment investigates the behaviour of the considered
$p_{dt}$ predictors for different time instants, $T$, and autocorrelation
parameters, $\rho$. Specifically, it calculates BP-plug-in, BP, plug-in and
EBP. Given the computational burden presented by the BPs (and EBPs) of the
target parameter and of the two random effects under Model ST1, the simulation
considers their approximated versions (see Section 3). The first predictor,
BP-plug-in, is obtained from (3.7) by using the theoretical vector of model
parameters $\theta$. The BPs and EBPs are approximated by generating $S1=500$
random variables $\hbox{\boldmath$v$}_{1}^{(s_{1})}$ and $S2=700$ random
variables $v_{2,\ell\tau}^{(s2)}$. For the empirical predictors (plug-in and
EBP), the model parameters are estimated by using the second option in MM,
since it has presented better results in the previous simulation experiment.

Table B.3 presents the average across domains and time instants of the biases
and the RMSEs (both $\times 10^{2}$) for BP-plug-in, BP, plug-in and EBP. BP
and EBP are more competitive than the respective plug-in. Specially in bias,
where a significant improvement is achieved. The obtained results also suggest
that $T$ does not affect too much the results and that variance is the most
important term of the MSE since bias is much smaller than RMSE.

Table B.3: Bias (B) and root mean squared error (RMSE) of the BP-plug-in, BP, plug-in and EBP of $p_{dt}$ under Model ST1 (both $\times 10^{2}$). |  | $T=4$ |  | $T=8$
---|---|---|---|---
$\rho$ | Predictors | B | RMSE |  | B | RMSE
0.1 | BP-plug-in | 0.3350 | 2.8325 |  | 0.3318 | 2.7739
| BP | 0.0701 | 2.7581 |  | 0.0689 | 2.7165
| plug-in | 0.3251 | 2.8248 |  | 0.3261 | 2.7710
| EBP | 0.0699 | 2.7645 |  | 0.0687 | 2.7203
0.3 | BP-plug-in | 0.3423 | 2.8445 |  | 0.3217 | 2.7840
| BP | 0.0711 | 2.7655 |  | 0.0674 | 2.7306
| plug-in | 0.3406 | 2.8354 |  | 0.3194 | 2.7851
| EBP | 0.0717 | 2.7723 |  | 0.0677 | 2.7342
0.5 | BP-plug-in | 0.3363 | 2.8939 |  | 0.3311 | 2.8125
| BP | 0.0692 | 2.8129 |  | 0.0829 | 2.7562
| plug-in | 0.3290 | 2.8830 |  | 0.3282 | 2.8169
| EBP | 0.0699 | 2.8200 |  | 0.0824 | 2.7595

The system of MM nonlinear equations (2.3) is solved by using the R nleqslv
package. The mvtnorm package is also used to generate samples of a SAR(1)
process. The computational burden of the first option in MM is much higher.
Taking $T=4$, the average runtime of the first option is $60.4$ seconds, while
for the second option is $0.1$ seconds. On the other hand, regarding the
computational burden of the $p_{dt}$ predictors, the EBP is faster than the
plug-in. The reason is because the proposed plug-in predictor requires the
calculation of two EBPs ($\hat{\hbox{\boldmath$v$}}_{1}$ and
$\hat{\hbox{\boldmath$v$}}_{2}$). The average runtimes are $210.7$ seconds for
the EBP and $320.5$ seconds for the plug-in.

## References

Baldermann, C., Salvati, N., and Schmid, T. (2016). Robust small area
estimation under spatial non-stationarity. Discussion Paper, School of
Business and Economics: Economics, 22, N. 2016/5.

Boubeta, M., Lombardía, M. J., W., Marey-Pérez, M., Morales, D. (2015).
Prediction of forest fires occurrences with area-level Poisson mixed models.
Journal of Environmental Management, 154, 151-158.

Boubeta, M., Lombardía, M. J., and Morales, D. (2016). Empirical best
prediction under area-level Poisson mixed models. Test, 25, 548-569.

Boubeta, M., Lombardía, M. J., and Morales, D.(2017). Poisson mixed models for
studying the poverty in small areas. Computational Statistics and Data
Analysis, 107, 32-47.

Boubeta, M., Lombardía, M. J., W., Marey-Pérez, M., Morales, D. (2019).
Poisson mixed models for predicting number of fires. International Journal of
Wildland Fire, 28, 3, 237-253.

Chandra, H., Salvati, N., Chambers, R., and Tzavidis, N. (2012). Small area
estimation under spatial nonstationarity. Computational Statistics and Data
Analysis, 56, 2875-2888.

Chandra, H., Salvati, N., and Chambers, R. (2015). A spatially nonstationary
Fay-Herriot model for small area estimation. Journal of Survey Statistics and
Methodology, 3, 109-135.

Chandra, H., Salvati, N., and Chambers, R. (2017). Small area prediction of
counts under a non-stationary spatial model. Spatial Statistics, 20, 30-56.

Chandra, H., Salvati, N., and Chambers, R. (2018). Small area estimation under
a spatially non-linear model. Computational Statistics and Data Analysis, 126,
19-38.

Choi, J., Lawson, A.B., Cai, B., Hossain, M.M. (2011). Evaluation of Bayesian
spatiotemporal latent models in small area health data. Environmetrics, 22, 8,
1008-1022.

Cressie, N. (1993). Statistics for spatial data. Wiley, New York.

Esteban, M. D., Morales, D., Pérez, A., L. Santamaría (2012). Small area
estimation of poverty proportions under area-level time models. Computational
Statistics and Data Analysis, 56, 10, 2840-2855.

Esteban, M. D., Morales, D., Pérez, A. (2016). Area-level spatio-temporal
small area estimation models. In Analysis of poverty data by small area
estimation (ed M. Pratesi). John Wiley and Sons, Ltd, Chichester, UK.

Esteban, M. D., Lombardía, M.J., López-Vizcaíno, E., Morales, D., Pérez A.
(2020). Small area estimation of proportions under area-level compositional
mixed models. TEST. DOI: 10.1007/s11749-019-00688-w.

González-Manteiga, W., Lombardía, M.J., Molina, I., Morales, D. and
Santamaría, L. (2008). Analytic and bootstrap approximations of prediction
errors under a multivariate Fay-Herriot model. Computational Statistics and
Data Analysis, 52, 5242-5252.

González-Manteiga, W., Lombardía, M. J., Molina, I., Morales, D., Santamaría,
L. (2010). Small area estimation under Fay-Herriot models with nonparametric
estimation of heteroscedasticity. Statistical Modelling, 10, 215-239.

Jiang, J. (1998). Consistent estimators in generalized linear models. Journal
of the American Statistical Association, 93, 720-729.

López-Vizcaíno, E., Lombardía, M.J., Morales, D. (2013). Multinomial-based
small area estimation of labour force indicators. Statistical Modelling, 13,
153-178.

López-Vizcaíno, E., Lombardía, M.J. and Morales, D. (2015). Small area
estimation of labour force indicators under a multinomial model with
correlated time and area effects. Journal of the Royal Statistical
Association, series A, 178, 535-565.

Marhuenda, Y., Molina, I., and Morales, D. (2013). Small area estimation with
spatio-temporal Fay-Herriot models. Computational Statistics and Data
Analysis, 58, 308 - 325. The Third Special Issue on Statistical Signal
Extraction and Filtering.

Meteogalicia. Xunta de Galicia. In web site:
https://www.meteogalicia.gal/observacion/informesclima/informesIndex.action.

Molina, I., Salvati, N., and Pratesi, M. (2009). Bootstrap for estimating the
MSE of the spatial EBLUP. Computational Statistics, 24, 441-458.

Moura, F. A. S. and Migon, H. S. (2002). Bayesian spatial models for small
area estimation of proportions. Statistical Modelling, 2(3), 183-201.

Opsomer, J. D., Claeskens, G., Ranalli, M. G., Kauermann, G., and Breidt, F.
J. (2008). Nonparametric small area estimation using penalized spline
regression. Journal of the Royal Statistical Society, Series B, 70, 265-286.

Pereira, L.N., Coelho, P.S. (2012). Small area estimation using a spatio-
temporal linear mixed model. REVSTAT - Statistical Journal, 10, 285-308.

Petrucci, A. and Salvati, N. (2006). Small area estimation for spatial
correlation in watershed erosion assessment. Journal of Agricultural,
Biological, and Environmental Statistics, 11, 169-172.

Pratesi, M. and Salvati, N. (2008). Small area estimation: the EBLUP estimator
based on spatially correlated random area effects. Statistical Methods and
Applications, 17, 113-171.

Särndal, C., Swensson, B., Wretman. J. (1992). Model assisted survey sampling.
Springer

Singh, B., Shukla, G., and Kundu, D. (2005). Spatio-temporal models in small
area estimation. Survey Methodology, 31, 183-195.

Sugasawa, S., Kawakubo, Y., and Ogasawara, K. (2015). Geographically weighted
empirical Bayes estimation via natural exponential family. Discussion Paper
No. 2015-01. Tokyo Institute of Technology.

Ugarte, M. D., Ibáñez, B., and Militino, A. F. (2006). Modelling risks in
disease mapping. Statistical methods in medical research, 15, 21-35.

Ugarte, M. D., Goicoa, T., and Militino, A. F. (2010). Spatio-temporal
modeling of mortality risk using penalized splines. Environmetrics, 21,
270-289.

You, Y. and Zhou, Q. M. (2011). Hierarchical Bayes small area estimation under
a spatial model with application to health survey data. Survey Methodology,
37, 25-37.