Tingting Huang

Gilbert Saporta

Huiwen Wang

Shanshan Wang
email: sswang@buaa.edu.cn

A Robust Spatial Autoregressive Scalar-on-Function Regression with t-distribution

Keywords: EM algorithm, FPCA, Functional linear model, Spatial autoregressive model, Spatial dependence, t-distribution

Modelling functional data in the presence of spatial dependence is of great practical importance as exemplified by applications in the fields of demography, economy and geography, and has received much attention recently. However, for the classical scalar-on-function regression (SoFR) with functional covariates and scalar responses, only a relatively few literature is dedicated to this relevant area, which merits further research. We propose a robust spatial autoregressive scalar-on-function regression (RSSoFR) by incorporating a spatial autoregressive parameter and a spatial weight matrix into the SoFR to accommodate spatial dependencies among individuals. The t-distribution assumption for the error terms makes our model more robust than the classical spatial autoregressive models under normal distributions. We estimate the model by firstly projecting the functional predictor onto a functional space spanned by an orthonormal functional basis and then presenting an expectation-maximization (EM) algorithm. Simulation studies show that our estimators are efficient, and are superior in the scenario with spatial correlation and heavy tailed error terms. A real weather dataset demonstrates the superiority of our model to the SoFR in the case of spatial dependence.

INTRODUCTION

Functional data are high-dimensional structured data that vary over a continuous domain. Examples are datasets recorded densely over time, space or time-space, like weather data, stock market data, trajectory data, diffusion tensor imaging (DTI) data and mass spectrometry data. Functional data analysis (FDA), as a new area of statistics, develops statistical methods to analyze information within functional data and has been applied in many subject areas, such as biology, medical sciences, meteorology, econometrics, finance, chemometrics and geophysics (Ramsay andSilverman (2002, 2005)). One of the most important tools in FDA is functional regression, including scalar-on-function regression (SoFR), function-on-scalar regression and function-onfunction regression. In particular, the first one, i.e., the classical SoFR [START_REF] Ramsay | Some tools for functional data analysis[END_REF]; [START_REF] Hastie | a statistical view of some chemometrics regression tools[END_REF]; [START_REF] Hall | Methodology and convergence rates for functional linear regression[END_REF]) with scalar responses and functional covariates, is of great interest and has been extensively studied.

Let Y be a centered scalar response variable, and X(t) be a second-order stochastic process on a compact interval Γ, E(X(t)) = 0 and E( Γ X 2 (t)dt) < ∞. A SoFR associates the functional predictor X(t) to the scalar response Y by where β(t) is the unknown slope function, and is the random error term, which is independent of X(t) and has zero mean and finite variance.

Regarding model (1.1), there has been a great deal of literature dedicated to its extensions, see [START_REF] Morris | Functional regression[END_REF] and [START_REF] Reiss | Methods for scalar-on-function regression[END_REF] for reviews. According to the roles the functional predictor X(t) plays in regression, the variants of the SoFR can be categorized into three groups, linear functional predictor regression, nonlinear functional predictor regression and nonparametric functional regression. The first group involves the generalized SoFR [START_REF] Marx | Generalized linear regression on sampled signals and curves: A P-spline approach[END_REF]; [START_REF] James | Generalized linear models with functional predictors[END_REF]), the multi-level SoFR [START_REF] Crainiceanu | Generalized Multilevel Functional Regression[END_REF]; [START_REF] Goldsmith | Longitudinal penalized functional regression for cognitive outcomes on neuronal tract measurements[END_REF]), the functional mixture regression [START_REF] Fang | Functional mixture regression[END_REF]) and the partial SoFR [START_REF] Shin | Partial functional linear regression[END_REF]). For the second team, multivariate nonlinear regressions with scalar regressors have been extended to the functional predictor case, such as the functional quadratic regression [START_REF] Yao | Functional quadratic regression[END_REF]), the single-functional index model [START_REF] Ait-Saïdi | Cross-validated estimations in the single-functional index model[END_REF]), the multiple-index model with functional covariates [START_REF] James | Functional adaptive model estimation[END_REF]) and the continuously additive model [START_REF] Müller | Continuously additive models for nonlinear functional regression[END_REF]). The last group, the nonparametric paradigm of the SoFR, is initially studied by [START_REF] Ferraty | Nonparametric functional data analysis: theory and practice[END_REF]. A key assumption in the aforementioned literature is that individuals are mutually independent, which may crumble. Yet researches concerning the SoFR with spatial dependence are relatively rare.

In reality, however, spatially correlated data are very common. For example, unemployment data [START_REF] Topa | Social Interactions, Local Spillovers and Unemployment[END_REF]) and housing prices data show spatial correlation, annual precipitation of a city relates to it's neighboring cities'. The characteristic of these data is that the dependency of two spatial units is determined by the distance between them.

This distance can be either Euclidean or more general [START_REF] Isard | General theory: social, political, economic and regional. Cambridge[END_REF]), for example, social distance, policy distance and economic distance. In spatial statistics and spatial economics, an abundance of publications are devoted to analyze data with spatial dependence [START_REF] Cressie | Statistics for spatio-temporal data[END_REF]; [START_REF] Schabenberger | Statistical methods for spatial data analysis[END_REF]; [START_REF] Anselin | Spatial Econometrics: Methods and Models[END_REF]; [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF]; [START_REF] Wang | Spatial partial least squares autoregression: Algorithm and applications[END_REF]). Discrete entities (areal data) and continuous surfaces (point-referenced data) are the two main research objects (a distinction between them refers to [START_REF] Anselin | Under the hood issues in the specification and interpretation of spatial regression models[END_REF]). For the last few years, there are many works fitting spatially correlated functional data, to name a few, [START_REF] Zhang | Functional concurrent linear regression model for spatial images[END_REF], [START_REF] Menafoglio | Statistical analysis of complex and spatially dependent data: a review of object oriented spatial statistics[END_REF], [START_REF] Nerini | Cokriging for spatial functional data[END_REF], [START_REF] Zhang | Functional CAR models for large spatially correlated functional datasets[END_REF], [START_REF] Giraldo | Spatial prediction of a scalar variable based on data of a functional random field[END_REF] and [START_REF] Aguilera-Morillo | Prediction of functional data with spatial dependence: a penalized approach[END_REF]. We find these papers are either interested in point-referenced data (for which kriging methods are commonly used) or related to regressions with functional responses. Articles concerning fitting functional regressors on scalar responses for discrete entities are relatively scanty. In this study, our interest lies in the SoFR for areal data.

To better illustrate our motivation for implementing this research and its significance, here is an application instance. We collected weather data from the China Meteorolog-ical Yearbook covering the period between 2005 and 2007. These data record monthly mean temperatures and monthly total precipitation in 34 major cities in China. Our aim is to investigate the effect of temperature on precipitation over these three years.

The SoFR (1.1) is popularly used to study this problem (see [START_REF] Ramsay | Functional Data Analysis[END_REF]). In the data preprocessing, we average the monthly temperatures and the monthly precipitation over 3 years to get their means, and add up the mean monthly total precipitation over 12 months to obtain the mean annual rainfall. The scalar response is the logarithm of the mean annual total precipitation, and the functional covariate is the mean monthly temperature. After applying the SoFR directly to the weather data, we want to examine whether there is spatial dependence among the residuals. The Moran's I test statistic [START_REF] Cliff | Testing for spatial autocorrelation among regression residuals[END_REF]) is employed (in Section 5

we explained why the weather dataset can be regarded as areal data instead of pointreferenced data). We find the resulting value of the Moran's I statistic is 0.48, and the p-value is smaller than 0.001, which indicates significant correlation exists among the residuals of the SoFR.

Figure 1: The Moran Scatterplot of the residuals of the SoFR.

For further illustration, we display the Moran Scatterplot of the residuals in Figure 1, from which it can be seen that there is approximately linear relationship between the residuals and the spatially lagged residuals (i.e., and nW , which are defined the same as those in (2.1)) of the SoFR, which means the SoFR may be not appropriate for spatially dependent data. This motivates the incorporation of spatial correlation in the analyses of the SoFR. A detailed analysis of this weather dataset and the construction of the weight matrix W can be found in Section 5.

When the predictor is scalar, instead of functional, linear regression models including spatial autoregressive (SAR) model, spatial error model (SEM) and spatial Durbin model (SDM) are frequently used to accommodate spatial dependence for areal data [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF]). It is straightforward to borrow concepts from the spatial models for modelling functional data with lattice structure. As the SAR model is representative among these spatial models and has been widely studied, we mainly focus on the SAR model and formulate our new model. Applications of the SAR model can be seen in [START_REF] Case | Spatial patterns in household demand[END_REF], [START_REF] Topa | Social Interactions, Local Spillovers and Unemployment[END_REF], and [START_REF] Olubusoye | Modelling road traffic crashes using spatial autoregressive model with additional endogenous variable[END_REF], among others. Estimation methods for the SAR model are described in [START_REF] Ord | Estimation methods for models of spatial interaction[END_REF], [START_REF] Lee | Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models[END_REF], [START_REF] Kelejian | A Generalized Moments Estimator for the Autoregressive Parameter in a Spatial Model[END_REF], [START_REF] Lee | GMM and 2SLS estimation of mixed regressive, spatial autoregressive models[END_REF] and [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF]. In an SAR model, a spatial weight matrix is employed to quantify adjacent relations among the observations, and an unknown spatial autoregressive parameter ρ is used to reflect the strength of spatial dependence. We take the advantages of the SAR model and incorporate the spatial correlation into the SoFR model using a spatial autocorrelation parameter and a weight matrix.

We noticed that in the above spatial models, the residual terms are generally presumed to follow normal distributions. This assumption is not realistic, as in applied problems we often expect tails of the error distribution are longer than the those of the normal distributions. Besides, when using the ordinary least squares (OLS) or maximum likelihood estimation method to estimate the spatial models, the estimated parameters can be largely affected by atypical points in dataset. Thus there is a need to consider non-normality error term in the spatial linear models. Fortunately, t-distribution, which has thick tails, provides an alternative. So we assume the residuals follow tdistributions in our model. The proposed new model is named the robust spatial autoregressive scalar-on-function regression (RSSoFR).

In this article, we presented a robust spatial functional linear model (RSSoFR), which extends the independent SoFR to the spatial scenario where the responses are spatially correlated. And by supposing the error terms follow t-distributions, the new model is more robust than the classical spatial models under normal distribution assumptions.

To estimate the RSSoFR, we firstly project the functional predictor onto a smaller functional space spanned by orthonormal basis functions, then propose an expectationmaximization (EM) algorithm to handle the spatial parameters. Two functional bases, the functional partial-least squares (FPLS) basis and the functional principal component (FPC) basis are considered. We also use a spatial cross-validation approach to select models. The article is organized as follows. In Section 2, we formulate the new model. The proposed estimation method is constructed based on orthonormal functional basis and the EM algorithm in Section 3. The finite-sample performances of the proposed estimators are evaluated through simulation studies in Section 4. Finally, in Section 5, we use a real dataset to document the usefulness of this methodology. We conclude the article with a discussion in Section 6.

MODEL SPECIFICATION

2.1 Spatial autoregressive (SAR) model [START_REF] Ord | Estimation methods for models of spatial interaction[END_REF] proposed a spatial autoregressive (SAR) model with parsimonious parameters. Assume {(x i , y i )} n i=1 are observed from n spatial units on a lattice, and denote

x = x 1 , x 2 , • • • , x n , y = (y 1 , y 2 , • • • , y n ) . The SAR model is, y = ρW y + xβ + , i ∼ N (0, σ 2 ) (2.1)
where W = (w ii ) n×n is a pre-specified spatial weight matrix, in which w ii represents the weight between units i and i , and = ( 1 , • • • , n ) is the noise term, which is independent of x, independently and identically follow normal distributions N (0, σ 2 ).

ρ and β are the parameters to be estimated. Here, the scalar parameter ρ ∈ (-1, 1).

In model (2.1), the spatial weight matrix W is exogenous, with each entry w ii assigned a value according to the contiguity or distances between units i and i in different contexts. For spatial contiguity matrices, the value of w ii is binary,

w ii =     
1, i and i are neighboring 0, i and i are not neighboring .

The assessment of neighbors relies on a known map indicating spatial arrangement of points, which can be regular or irregular. On a regular grid, units are neighbors if they share a border (rook case), a vertex (bishop case), or share either a border or a vertex (queen case) (for more details refer to [START_REF] Anselin | Spatial Econometrics: Methods and Models[END_REF]). In an irregularly spaced case, units are also neighboring when they have common edges. In the case of knowing geographical locations of spatial units, w ii can be functions of distance d ii between two points. Popular choices are the inverse distance (IV) and the negative exponential (NE) of distance,

w IV ii = 1 d ii , w N E ii = e -d ii . (2.2)
When the units are areal data, boundary length information can be also used to form w ii , see for example [START_REF] Decey | A review of measure of contuity for two and k-color maps[END_REF].

We mention some more general spatial matrices. In a social network, w ii is normally set to 1 if persons i and i are friends, 0 otherwise. Note that W is not necessarily a symmetric matrix. For example on a social media Sina Weibo, person i is a follower of person i but i does not follow i, we have w ii = 1, w i i = 0. Some economic factors such as the GDP (Gross Domestic Product) and income can be also used to establish W , for example, [START_REF] Case | Budget spillovers and fiscal policy interdependence: Evidence from the states[END_REF] used

w ii = 1 |IN C i -IN C i | ,
where IN C i is per capita income in state i. More information for W refers to [START_REF] Isard | General theory: social, political, economic and regional. Cambridge[END_REF] andAnselin (1998).

In general, we standardize the spatial matrix W to be a row-normalized matrix; in this matrix, the summation of the row elements is unity, and the entries on the diagonal are zeros. About the interpretations of model (2.1), refer to [START_REF] Anselin | Under the hood issues in the specification and interpretation of spatial regression models[END_REF].

Robust spatial autoregressive scalar-on-function regression (RSSoFR)

Following [START_REF] Qu | Estimating a spatial autoregressive model with an endogenous spatial weight matrix[END_REF], consider the spatial processes located on an unevenly spaced lattice D ⊆ R d , d ≥ 1. And we observe {(x i (t), y i )} n i=1 from n spatial units on D. Here, x i (t)s are square integrable second-order stochastic processes defined on a compact set Γ. Without loss of generality, we presume Γ is the unit,

i.e. t ∈ [0, 1]. Denote x(t) = x 1 (t), x 2 (t), • • • , x n (t) , we formulate the RSSoFR as y = ρW y + 1 0 x(t)β(t)dt + , i ∼ t(ν) (2.3) where W is pre-defined as the W in model (2.1), = ( 1 , • • • , n ) independently and
identically follow t-distributions with freedom ν. Here, ν, β(t) and ρ ∈ [0, 1) are the parameters to be estimated.

As mentioned previously, ρ is a scale parameter that reflects strength of the impacts from neighbours. Greater values of ρ indicate that y i is more strongly affected by its neighbours. The W matrix can be pre-formed as the W pre-specified in model (2.1).

We use the inverse distance to construct W in Section 5 for real data analysis. The proposed RSSoFR is more general than the SoFR and the SAR model.

• When ρ = 0, the RSSoFR reduces to an SoFR.

• When x(t) is free of t, the RSSoFR degenerates into an SAR model with tdistributions.

To provide better insight onto the new model, we reformulate equation ( 2.3) as the following equivalent expression,

y = (I n -ρW ) -1 1 0 x(t)β(t)dt + (I n -ρW ) -1 ,
which shows how y is generated. We can know the mean of y given

x(t) is E y|x(t) = (I n -ρW ) -1 1 0 x(t)β(t)dt.
Thus the spatial process of y is not stationary, i.e., the mean of y depends on the spatial units' locations through the functional covariate x(t).

The error terms (I n -ρW ) -1 shows the residuals of y i s are spatially correlated. As each i follows a t-distribution, the MLE method can not be used. We put forward an EM algorithm to estimate the parameters in (2.3).

ESTIMATION METHOD

In this section, we first expand the functional predictor by orthogonal functional basis, then present an EM algorithm to obtain the estimators of the spatial autocorrelation parameter ρ, the slope function β(t), and the freedom parameter ν in the RSSoFR (2.3).

Expand the functional predictor by orthonormal basis

Note that, before estimating the parameters of our model, data representation methods, such as smoothing and interpolation, should be used to convert discretely recorded data x i (t j ) to curves x i (t). And x(t), y are centered in advance.

A standard method to handle the functional predictor in the SoFR is expressing the curves by a linear combination of orthonormal basis functions. We also firstly expand the functions in (2.3) so that the functional term 1 0 x(t)β(t)dt can be represented by finite real vectors. Suppose {φ j (t)} ∞ j=1 is a standard orthogonal basis of L 2 space, which is composed of square integrable functions. Then x i (t) and β(t) can be written

as x i (t) = ∞ j=1 a ij φ j (t) and β(t) = ∞ j=1 b j φ j (t), where a ij = 1 0 x i (t)φ j (t)dt, b j = 1 0 β(t)φ j (t)dt. Moreover, 1 0 x i (t)β(t)dt = ∞ j=1 a ij b j by orthogonality. Denote a j = (a 1j , a 2j , . . . , a nj ) , we can rewrite model (2.3) as y = ρW y + ∞ j=1 a j b j + , (3.1)
which is easier to cope with under the numerical vectors framework. In practice, the estimation of β(t) is an ill-posed problem. Regularization procedure is thus needed.

Some authors add a penalty term in the objective function to put constraints on β(t) [START_REF] Cardot | Spline estimators for the functional linear model[END_REF]; [START_REF] Crambes | Smoothing splines estimators for functional linear regression[END_REF]). Here, we project x(t) and β(t) onto a finite dimensional space spanned by m basis functions {φ j (t)} m j=1 . Therefore, expression (3.1) has the truncated form

y ≈ ρW y + m j=1 a j b j + . (3.2)
Regarding the orthonormal functional basis φ j (t)s, we can use the Fourier basis, the 2012)). Here, we both introduce the FPC basis and the FPLS basis in the process of estimation. As these two bases are adaptive to properties of data, whereas other bases independent of data will make it hard to decide which m terms should be included in the truncated RSSoFR (3.2).

Functional principal component (FPC) basis

Let K(s, t) denote the covariance function of X(t), i.e. K(s, t) = Cov(X(t), X(s)). By Mercer's theorem, the spectral decomposition of K(s, t) is then K(s, t) = ∞ j=1 λ j ϕ j (s) ϕ j (t), where λ 1 > λ 2 > • • • > 0 are eigenvalues and {ϕ j (t)} ∞ j=1 are the corresponding orthogonal eigenfunctions. According to the Karhunen-Loève expansion, X(t) can be expanded as

X(t) = ∞ j=1 a j ϕ j (t)
, where a j s are uncorrelated random variables with mean zero and variance E(a 2 j ) = λ j , a j = 1 0 X(t)ϕ j (t)dt.

For the observations

{y i , x i (t)} n i=1 , the empirical version of K(s, t) is K(s, t) = 1 n n i=1 x i (s) x i (t)
. Moreover, it can be shown that K(s, t) = n j=1 λj φj (s) φj (t), where λj and φj (t) are the estimators of ϕ j (t) and λ j , respectively. For the ith observation x i (t), the estimator of a ij is then âij = 1 0 x i (t) φj (t), and x i (t) can be written as x i (t) = where âj = (â 1j , â2j , • • • , ânj ) .

Functional partial-least squares (FPLS) basis

In the analysis of numerical data, the partial-least squares (PLS) regression is an efficient alternative to the principal component regression (PCR) by taking the response variable y onto account. Like the PCR, the PLS regression also provides a set of orthogonal basis functions. De Jong (1993) has proved that with the same number of components, the PLS regression has a better fit than the PCR. Especially when a large proportion of the variations of x i (t)s do not explain the responses y, the PLS regression has competitive advantages over the FPC regression.

Here, to form the FPLS basis functions, we first neglect the autoregressive term ρW y of (2.3), which means building the basis without considering spatial correlation. Then, use an iterative process introduced by [START_REF] Preda | PLS regression on a stochastic process[END_REF] to get the basis.

Supposing we construct m FPLS basis functions, the main steps are,

(1) Begin from l = 1. Also write x l (t) = x(t), y l = y.

(2) Obtain a square integrable weight function ω l (t) evaluated by ω l (t) = E[y l x l (t)] E[y l x l (t)] , which maximizes Cov

ω l (t) =1 y l , 1 0 x l (t)ω l (t)dt , where • is the norm, i.e., ω l (t) = 1 0 ω 2 l (t)dt.
(3) Regress x l (t) and y l on a l , separately, where

a l = 1 0 x l (t)ω l (t)dt. That is x l (t) = p l (t)a l + x l (t) and y l = q l a l + y l , where p l (t) = E[x l (t)a l ] a l 2 , q l = E[y l a l ] a l 2 . (4) Stop when l = m. Otherwise, set x l+1 (t) = x l (t)
, y l+1 = y l , and back to step 2.

The set of weight functions {ω j (t)} m j=1 is the FPLS basis. And the orthogonal components a 1 , . . . , a m are the regressors a j of the truncated model (3.2). In practice, we compute E[y l x l (t)], E[x l (t)a l ] and E[y l a l ] by their empirical versions. Denote âj as the estimator of a j , we have the same expression as (3.3), which is the truncated empirical form of (3.1).

An EM algorithm for the truncated RSSoFR

In this subsection we focus on (3.3) and propose an EM algorithm to estimate ρ, b = ( b1 , b2 , • • • , bm ) and ν. Defining A = (â ij ) n×m , the expression (3.3) can be written as y ≈ ρW y + Ab + .

(3.4)

It is not easy to directly write out the probability density function of y based on the noise term as = ( 1 , • • • , n ) independently and identically follow t-distributions.

However, we know the fact that a t-distribution can be regarded as a scale mixture of normal distributions, i.e., if u i follows a gamma distribution f

(u i ) = 1 Γ( ν 2 ) ( ν 2 ) ν 2 (u i ) ν 2 -1 e -ν 2 u i ,
and i | u i ∼ N (0, σ 2 u i ), then the marginal density f ( i ) of i is a t-distribution with
freedom ν and scale parameter σ,

f ( i ) = Γ( ν+1 2 )σ -1 (πν) 1 2 Γ( ν 2 ){1 + 2 σ 2 ν } 1 2 (ν+1) .
Thus by introducing latent variables u = (u 1 , u 2 , . . . , u n ) , which are independent of x(t) and independently and identically distributed (i.i.d.), we put forward an EM algorithm to estimate (3.4).

Given u, (3.4) can be expressed as

y | u ≈ (I n -ρW ) -1 (Ab + σU -1 2 η) η ∼ N (0, I n ), U = diag(u)
where η = (η 1 , η 2 , . . . , η n ) . Then the complete log-likelihood function for y and u is ln

L(θ; y, u) = ln f (y|u) + n i=1 ln f (u i ) = - n 2 ln(2π) - 1 2 ln |Σ y | - 1 2 (y -µ y ) Σ -1 y (y -µ y ) + n i=1 ln f (u i ), (3.5)
where θ = (ρ, b, σ 2 , ν), µ y and Σ y are the mean and the covariance of y, respectively,

µ y = (I n -ρW ) -1 Ab, Σ y = (I n -ρW ) -1 • σ 2 • U -1 • (I n -ρW ) -1 .
Denote the constant term in (3.5) by C, the complete log-likelihood function is

ln L(θ; y, u) =C + 1 2 n i=1 ln u i - n 2 ln σ 2 + ln |I n -ρW | - U 2σ 2 + n i=1 ( ν 2 -1) ln u i - ν 2 u i -ln Γ( ν 2 ) + ν 2 ln ν 2 .
(3.6) Where = y -ρW y -Ab. Observing the above, if we known U , the estimators of ρ, b, σ 2 are directly obtained by the quasi-maximum likelihood estimation method introduced by [START_REF] Lee | Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models[END_REF]. And the freedom ν can be also estimated by maximizing the last term in (3.6), i.e.

n i=1 ( ν 2 -1) ln u i - ν 2 u i -ln Γ( ν 2 ) + ν 2 ln ν 2 .
However, u is missing data, and next we present the EM algorithm.

In E-step, given the kth estimate of θ, the computation of E ln L(θ; y, u|y, θ (k) )

comes down to compute E(u|y, θ (k) ) and E( n i=1 ln u i |y, θ (k) ). Note that the fact that the conditional distribution of u i given i is also a gamma distribution G(ν 1 , ν 2 ),

where

ν 1 = 1 2 (ν + 1), ν 2 = 1 2 (ν + 2 i σ 2
), and we have

E(u i | i ) = ν+1 ν+( 2 i )/(σ 2 ) . Hence E(u|y, θ (k) ) = ν (k) + 1 ν (k) + δ (k) 1 , ν (k) + 1 ν (k) + δ (k) 2 , . . . , ν (k) + 1 ν (k) + δ (k) n , δ (k) i = ( 2 i ) (k) (σ 2 ) (k) .
We also have

E n i=1 ln u i |y, θ (k) = n i=1 ψ ν (k) + 1 2 -ln ν (k) + δ (k) i 2 , where ψ(s) = ∂ ln Γ(s)
∂s is the Digamma function.

And in M-step, the (k + 1)th estimate of θ maximizes

E ln L(θ; y, u|y, θ (k) ) ∝ ν 2 E n i=1 ln u i |y, θ (k) - ν 2 τ n E(u|y, θ (k) ) -n ln Γ( ν 2 )+ n ν 2 ln ν 2 - n 2 ln(σ 2 ) + ln |I n -ρW | - E(U |y, θ (k) ) 2σ 2 ,
where τ n is an n-dimensional vector of ones, k) ) . Thus we can get ν (k+1) by maximizing

E(U |y, θ (k) ) = diag E(u|y, θ (
ν 2 E n i=1 ln u i |y, θ (k) - ν 2 τ n E(u|y, θ (k) ) -n ln Γ( ν 2 ) + n ν 2 ln ν 2 .
And obtain {ρ (k+1) , b (k+1) , (σ 2 ) (k+1) } by maximizing

n 2 ln(σ 2 ) + ln |I n -ρW | - E(U |y, θ (k) ) 2σ 2 . (3.7)
Observing the above if ρ is known, b (k+1) and (σ 2 ) (k+1) are respectively,

b (k+1) (ρ) = A E(U |y, θ (k) )A -1 A E(U |y, θ (k) )y -ρA E(U |y, θ (k) )W y (σ 2 ) (k+1) (ρ) = 1 n y -ρW y -Ab (k+1) (ρ) • E(U |y, θ (k) ) • y -ρW y -Ab (k+1) (ρ) .
(3.8)

Thus similar to the quasi-maximum likelihood estimation method [START_REF] Lee | Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models[END_REF]), we substitute b (k+1) (ρ) and (σ 2 ) (k+1) (ρ) onto (3.7) and drop the constant term. Then

ρ (k+1) is the solution of arg max ρ - n 2 ln (σ 2 ) (k+1) (ρ) + ln |I n -ρW | .
The optimization of the above refers to [START_REF] Ord | Estimation methods for models of spatial interaction[END_REF]. Once ρ (k+1) is obtained, replacing ρ with ρ (k+1) in (3.8) yields b (k+1) and (σ 2 ) (k+1) , respectively.

The estimator of θ is finally obtained by repeating E-step and M-step until max{|ρ (k+1) - (3.9)

ρ (k) |, |b (k+1) -b (k) |, |(σ 2 ) (k+1) -(σ 2 ) (k) |}
For clarity, we summarize the main steps of the estimation procedure in Algorithm 1.

Choosing the truncation parameter for the RSSoFR

We introduce two ways to determine the truncation parameter m. The first is the percentage of variance explained (PVE) for predictors, which is often adopted for the FPC basis. The second is based on cross validation, which is a more general approach.

Algorithm 1 Main steps of the estimation procedure 1: Represent the functional predictor and the slope function using the FPLS basis or the FPC basis. In this step, after an appropriate truncation parameter is given, the RSSoFR approximates an SAR model with t-distributions whose covariates are the component scores of x i (t)s, as shown in (3.4). 2: Determine the estimators for ρ, b, σ 2 , ν in (3.4) by the EM algorithm.

E-step : Calculate u (k+1) i and (ln u i ) (k+1) .

u (k+1) i = E(u i |y, θ (k) ) = ν (k) + 1 ν (k) + δ (k) i (ln u i ) (k+1) = E(ln u i |y, θ (k) ) = ψ ν (k) + 1 2 -ln ν (k) + δ (k) i 2 where ψ(s) is the Digamma function, δ (k) i = ( 2 i ) (k) (σ 2 ) (k) and (k) i is the ith element of y -ρ (k) W y -Ab (k) . M-step : Update the parameters {ν (k) , ρ (k) , b (k) , (σ 2 ) (k) }. ρ (k+1) = arg max ρ - n 2 ln (σ 2 ) (k+1) (ρ) + ln |I n -ρW | b (k+1) = A U (k+1) A -1 A U (k+1) y -ρ (k+1) A U (k+1) W y (σ 2 ) (k+1) = 1 n y -ρ (k+1) W y -Ab (k+1) • U (k+1) • y -ρ (k+1) W y -Ab (k+1) ν (k+1) = arg max ν n i=1 ν 2 (ln u i ) (k+1) - ν 2 u (k+1) i -ln Γ( ν 2 ) + ν 2 ln ν 2
where U (k+1) = diag (u

(k+1) 1 , • • • , u (k+1) n
) . 3: Reconstruct the estimator for β(t) in the RSSoFR. The slope function is constructed using the FPLS basis or the FPC basis mentioned in Step 1 and the coefficient b estimated in Step 2, as shown in (3.9). The other estimators are obtained directly from Step 2.

To compare the new RSSoFR with the SSoFR and the SoFR in the numerical experiments, we use the first method to determine m. If the PVE is set to be 80%, the truncation parameter m is subject to min l {( l j=1 λj )/( n j=1 λj ) ≥ 80%}.

In a non-spatial scenario, cross-validation method is commonly conducted by each time extracting an object from data as test set, and using the remaining as training set to predict the extracted one. In a spatial context, cross validation can be also applied but notice that the modification of spatial structure resulting from removing one unit should be considered. Denote the objects in training set as in-sample units, whose dependent variable Y s and independent variables X s are both observed, and the objects in test set as out-sample units, whose observed variable is explanatory variable X o and the response Y o is unknown, which need to be predicted. We call the problem of predicting Y o using X s , X o and Y s the out-of-sample prediction. [START_REF] Goulard | About predictions in spatial autoregressive models: optimal and almost optimal strategies[END_REF] discussed the out-of-sample prediction for the SAR model and concluded the "BP" predictor, which is based on Goldberger formula, behaves well among other existing predictors. The calculation formula of the BP predictor is

Ŷ BP o = Ŷ T C o -Q -1 o Q os (Y s -Ŷ T C s ), Ŷ T C = (I n -ρW ) -1 X β =   Ŷ T C s Ŷ T C o   ,
where Q is the precision matrix and

Q = 1 σ2 (I n -ρ(W + W ) + ρ2 W W ) =   Q s Q so Q os Q o   .
In all, we use a spatial cross-validation method, which employs the BP predictor, to select the number of components. Specifically, we extract one observation at a time from the dataset and use the remaining to predict it. When all the observations in the dataset have been extracted once, we define the following cross-validated prediction error (PE) for truncation parameter m,

P E(m) = 1 n n i=1 (y i -ŷm i ) 2 , (3.10)
where ŷm i is the predicted y i using m components and n -1 observations. For all the alternatives of m, we choose the one with the smallest PE.

SIMULATION STUDY

Several simulation studies are conducted to evaluate the finite-sample performances of the proposed estimators for ρ, ν, β(t). All of the computations were carried out in the R environment, and we used the R packages 'spdep', 'fda' and 'fda.usc' to implement the proposed procedure.

Because the spatial dependence and the t-distribution are of interest in this study, we compare the proposed RSSoFR with the following two models (a) The SoFR with i.i.d. t-distributions

y = 1 0 x(t)β(t)dt + , i ∼ t(ν). (4.1) (b) The SSoFR with i.i.d. normal distributions y = ρW y + 1 0 x(t)β(t)dt + , i ∼ N (0, σ 2 ). (4.2)
These three models are estimated all by firstly expanding the functional regressor using the FPC basis. Here we set the PVE to be identical to 80%. Then we have

(1) The SSoFR (4.2) approximates an SAR model with i.i.d. normal distributions, which can be estimated by the quasi-maximum likelihood estimation method [START_REF] Lee | Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models[END_REF]).

(2) The SoFR (4.1) approximates an ordinary multivariate linear model with i.i.d. tdistributions, which can be estimated by the EM algorithm [START_REF] Peel | Robust mixture modelling using the t distribution[END_REF]).

Regarding β(t) for the three models, they can be reconstructed by the basis functions as that showed by (3.9). Different degrees of spatial effects are considered, ρ = {0, 0.5, 0.8}. Note that when ρ = 0, the SSoFR (2.3) reduces to the SoFR (4.1).

As for the spatial scenario, we adopt the rook matrix by randomly apportioning n agents on a regular square grid of cells; each agent is located on a cell. In this context, if the grid has R rows and T columns, then the sample size n = R × T . Units that share an edge are neighbours. This definition ensures the units in the inner field of the grid have four neighbours, the units in the corners have two neighbours, and the units along the boarders have three neighbours. Therefore, the spatial matrix is an adjacent matrix with each entry w ii = 1 if units i and i are neighbours and w ii = 0 otherwise.

We set n = {10 × 15, 15 × 20, 20 × 30} in the simulation. The spatial weight matrix is row-normalized in all cases.

For the functional part of equation ( 2.3), we employ the same form as the functions in the SoFR by [START_REF] Hall | Methodology and convergence rates for functional linear regression[END_REF]. Specifically, we generate the simulation data

y = (y 1 , y 2 , • • • , y n ) using y = (I n -ρW ) -1 1 0 x(t)β(t)dt + 0.5 (1) i ∼ N (0, 1), i ∼ t(1), i ∼ t(3)
(2) ρ = 0, ρ = 0.5, ρ = 0.8

The functional predictor

x(t) = (x 1 (t), x 2 (t), • • • , x n (t)) is produced with values of x i (t) independently generated by x i (t) = 50 j=1 a j Z j ϕ j (t)
, where a j = (-1) j+1 j -γ/2 with γ = 1.1 and 2, respectively;

Z j ∼ U [- √ 3, √ 3
] and ϕ j (t) = √ 2 cos(jπt). Similarly, the coefficient function β(t) is generated according to β(t) = 50 j=1 b j ϕ j (t), where b 1 = 0.3 and b j = 4(-1) j+1 j -2 , j ≥ 2. Note that the eigenvalues of the covariance function K(u, v) play a vital role in determining the estimation accuracy of β(t). We consider two cases. In case 1, γ = 1.1, where the eigenvalues are well spaced and the slope function can be accurately estimated. In case 2, γ = 2, where the closely spaced eigenvalues can cause the estimator β(t) to display poor performance.

The experiment is repeated 500 times in each setting. And we suppose x i (t)s are observed on 101 equispaced points, i.e. t j = 0, 1 100 , • • • , 99 100 , 1. The selected number of components are {7, 8, 9} for γ = 1.1 and 2 for γ = 2. We assess the behaviour of the estimator ρ in terms of the mean bias and its standard deviation. Concerning β(t), we evaluate its performances in terms of the integrated mean-square error IM SE = 1 0 ( β(t) -β(t)) 2 dt. We summarize the estimation results for ρ in Table 1, β(t) in Table 2 andν in Table 3 better than that of the SSoFR. And when i follows N (0, 1) or t(3), ρ of the two models behave similarly.

(2) In Table 2 compare β(t) of the RSSoFR and the SoFR. In case ρ = 0, RSSoFR's β(t) approximates the SoFR's. But in case ρ = 0, β(t) of the RSSoFR has much smaller IMSE than that of the SoFR has. Moreover, as ρ increases, the differences in the IMSE between the RSSoFR and the SoFR also raise. Thus we conclude that overall our RSSoFR behaves better than the SoFR.

(3) In Table 2 observe β(t) of the RSSoFR and the SSoFR. We can find under normal distributions, our proposed method performs basically consistent with the SSoFR.

Whereas under t-distributions with thick tails, the RSSoFR produces better results for β(t) than the SSoFR. Thus the RSSoFR is more robust than the SSoFR.

(4) From Table 2, we can also find the IMSE of β(t) of the RSSoFR decreases with the increasing sample size. Besides, β(t) is more accurately estimated given γ = 1.1 than γ = 2 when the other simulation parameters are held equal, as the case in [START_REF] Hall | Methodology and convergence rates for functional linear regression[END_REF]. The estimation results for ν are presented in Table 3.

Figure 2: The estimated β(t) vs. the true β(t) when n = 600, γ = 1.1, ρ = 0, 0.5, 0.8 and i ∼ N (0, 1), t(1), t(3) respectively.

To show the performances of β(t) intuitively, for each case we select one result from 500 repetitions when n = 600, γ = 1.1, ρ = 0, 0.5, 0.8 and i ∼ N (0, 1), t(1), t(3). Figure 2 displays the estimated β(t) vs. the true value of β(t) for the RSSoFR, the SSoFR and the SoFR. From the figure, we can draw the similar conclusions for β(t) as those derived from Table 2. In addition, we have also implemented numerical experiments with the FPLS basis used for the functional regressor expansion, which leads to the where w ii is the weight between city i and i . We also build the SoFR as

y i = 1 0 x i (t)β(t)dt + i (5.2)
to enable comparison with the RSSoFR.

The spatial weight matrix W = (w ii ) n×n is formed by the reciprocal of the distance d ii between centers of two cities, i.e.,

w ii = 1 d ii ,
where d ii is computed using the haversine formula based on longitudes and latitudes of cities' centers. Besides, we considered two influence factors that may effect model (5.1)'s performances during the process of constructing W . The first is the threshold distance d 0 (in kilometers). We know if city i is very far from city i , the spatial dependence between them will be very small. Thus we set

w ii = 0 if d ii > d 0 .
The second is the number of the nearest neighbors k. For any city i, there exists competitions among weights w ii , i = 1, • • • , 34. Suppose city i has many close neighbors

i 1 , • • • , i n 0 , max{d ii 1 , • • • , d iin 0 } < d 0
and n 0 is relatively great, then it is very important to determine the value of k. As the smaller k is, the greater each w ii becomes, i.e., the stronger influence each neighbor i has on city i. 5. Figure 5 shows the related consequences. Firstly, we discuss the fitting results of the RSSoFR and the SoFR. Recall from Section 1 that significant spatial correlation presents among the residuals of the SoFR. We can see from Table 5 that the value of the Moran's I statistic of the residuals of the RSSoFR is small (-0.14), which suggests a majority of spatial correlation in the residuals has been removed. Figure 5 (top-left) shows the Moran Scatterplot of the residuals of the RSSoFR. As for the fitted values, the MSE of the RSSoFR is nearly half of that of the MSE of the SoFR, where M SE = 1 n n i=1 (ŷ i -y i ) 2 . Thus, the RSSoFR fits the mean total annual precipitation better than the SoFR. Figure 5 (bottomleft) presents boxplots of the residuals of the RSSoFR and the SoFR. Regarding the estimated parameters, ρ is 0.705, which is also very significant with its p-value smaller than 0.001. The slope functions of the RSSoFR and the SoFR are provided in Figure 5 (top-right). We can find the two estimated curves have similar shapes, and β(t) of the RSSoFR is smoother than that of the SoFR. We conclude that the precipitation is much more strongly influenced by the temperatures during spring and winter than in the other seasons. Then we focus on the prediction results. The predicted error of the RSSoFR is 0.076 while that of the SoFR is 0.114, which indicates the RSSoFR has better prediction performance than the SoFR. Figure 5 (bottom-right) presents boxplots of the prediction error under the RSSoFR and the SoFR.

. k = 1 k = 2 k = 3 k = 4 k = 5 k = 6 k = 7 k = 8 d 0 (

CONCLUSION AND DISCUSSION

The SoFR is popular in studies of links between a scalar response and functional predictors. However, the existing SoFR cannot address the dependencies in a crosssectional spatial scenario. We propose a robust spatial autoregressive scalar-on-function regression that incorporates a spatial lagged term into the SoFR to accommodate the spatial dependence and allows for thick tailed noise term. An estimation method based on basis expansion and EM algorithm is developed to obtain the estimators of the spatial autoregressive parameter and the slope function. Specifically, the FPC basis and the FPLS basis can be applied and the spatial cross-validation method is introduced to choose the truncation parameter. Our simulation study demonstrates the consistency of the proposed estimators. In particular, the new model performs better than the SoFR when the spatial correlation is present, and the SSoFR when the error tors via the bootstrap. There are some existing methods for estimating the standard error of an estimator obtained from the EM algorithm. For example, [START_REF] Louis | Finding the observed information matrix when using the em algorithm[END_REF] developed a formula for computing the observed information matrix in terms of the complete and missing information matrices, while the calculation of the missing information involves the conditional expectation of the outer product of the complete-data score vector; [START_REF] Jamshidian | Standard errors for em estimation[END_REF] proposed to use numerical differentiation to yield the Hessian matrix. After obtaining the standard errors, one may follow the analogous methods to perform the statistical inference for the proposed model. Investigations along this direction may be of interest and merit further research but are beyond the scope of the paper.

  FPC basis[START_REF] Dauxois | Asymptotic theory for the principal component analysis of a vector random function: Some applications to statistical inference[END_REF];[START_REF] James | Principal component models for sparse functional data[END_REF];[START_REF] Li | Uniform convergence rates for nonparametric regression and principal component analysis in functional/longitudinal data[END_REF]) and the FPLS basis[START_REF] Preda | PLS regression on a stochastic process[END_REF];Preda et al. (2007);Preda and Saporta (2007);[START_REF] Aguilera | Using basis expansions for estimating functional PLS regression: Applications with chemometric data[END_REF];Delaigle and Hall (

  n j=1 âij φj (t). Similarly, based on the estimated orthonormal functional basis { φj (t)} n j=1 , β(t) has the expansion β(t) = n j=1 bj φj (t) with bj = 1 0 β(t) φj (t)dt. Therefore, the sample counterpart of (3.2) is y ≈ ρW y + m j=1 âj bj + .(3.3)

  is within a given threshold. And β(t) can be reconstructed by the estimator of b, i.e. b = ( b1 , • • • , bm )

Figure 3 :

 3 Figure 3: The temperature curves of 34 cities in training set (left), the eigenvalues of the sample covariance function (right).

  Figure 4 presents the locations of 34 major cities on map of China based on the longitudes and the latitudes of the cities' centers. Because Urumchi and Lhasa are located far from the other cities (all the distances are greater than 1250 km), we remove their records from the weather dataset. Note that W is row-normalized after construction.

Figure 4 :

 4 Figure 4: The locations of 34 major cities on map of China based on the longitudes and the latitudes of the cities' centers.

Figure 3 (

 3 Figure 3 (right) shows the eigenvalues of the sample covariance function. The eigenvalues clearly decay quickly, and the first eigenvalue accounts for 98.6% of the total variance. Therefore, we consider using one component in the RSSoFR. As the FPLS basis behaves better than the FPC basis in this dataset, we employ the FPLS basis

Figure 5 :

 5 Figure 5: The Moran Scatter plot of the residuals of the RSSoFR (top-left), the estimated β(t) of the RSSoFR and the SoFR (top-right), the fitted error of the RSSoFR and the SoFR (bottom-left), and the predicted error of the RSSoFR and the SoFR (bottom-right), respectively.

Table 1 :

 1 . Examination of Table1-3 leads to the following conclusions. The biases and its standard deviations (in brackets) of ρ.

					γ = 1.1			γ = 2
	ρ	i	n	RSSoFR	SSoFR	SoFR	RSSoFR	SSoFR	SoFR
	0	normal 150	-0.0059	-0.0035	-	-0.0086	-0.0055	-
				(0.0887)	(0.0769)		(0.0969)	(0.0862)
			300	-0.0025	-0.0027	-	-0.0033	-0.0021	-
				(0.0608)	(0.0526)		(0.0646)	(0.0571)
			600	-0.0009	-0.0004	-	-0.0011	-0.0006	-
				(0.0429)	(0.0362)		(0.0479)	(0.0435)
	0	t(1)	150	0.0002	-0.0145	-	-0.0009	-0.0165	-
				(0.0266)	(0.0828)		(0.0291)	(0.0884)
			300	0.0004	-0.0092	-	-0.0003	-0.0090	-
				(0.0127)	(0.0600)		(0.0137)	(0.0552)
			600	-0.0006	-0.0042	-	-0.0003	-0.0009	-
				(0.0067)	(0.0345)		(0.0071)	(0.0338)
	0	t(3)	150	-0.0059	-0.0055	-	-0.0050	-0.0101	-
				(0.0801)	(0.0865)		(0.0862)	(0.0952)
			300	-0.0061	-0.0094	-	-0.0020	0.0000	-
				(0.0550)	(0.0596)		(0.0640)	(0.0708)
			600	-0.0021	-0.0011	-	-0.0026	-0.0033	-
				(0.0398)	(0.0467)		(0.0450)	(0.0511)
	0.5 normal 150	-0.0121	-0.0103	-	-0.0129	-0.0114	-
				(0.0733)	(0.0622)		(0.0775)	(0.0682)
			300	-0.0078	-0.0058	-	-0.0039	-0.0021	-
				(0.0500)	(0.0430)		(0.0554)	(0.0495)
			600	-0.0041	-0.0030	-	-0.0021	-0.0014	-
				(0.0362)	(0.0315)		(0.0404)	(0.0346)
	0.5	t(1)	150	-0.0016	-0.0192	-	0.0016	-0.0146	-
				(0.0220)	(0.0765)		(0.0217)	(0.0674)
			300	-0.0012	-0.0085	-	-0.0005	-0.0077	-
				(0.0111)	(0.0464)		(0.0096)	(0.0487)
			600	-0.0002	-0.0053	-	0.0000	-0.0034	-
				(0.0062)	(0.0364)		(0.0052)	(0.0296)
	0.5	t(3)	150	-0.0105	-0.0177	-	-0.0190	-0.0218	-
				(0.0712)	(0.0754)		(0.0785)	(0.0815)
			300	-0.0036	-0.0069	-	-0.0082	-0.0089	-
				(0.0467)	(0.0533)		(0.0504)	(0.0557)
			600	-0.0021	-0.0016	-	-0.0024	-0.0059	-
				(0.0309)	(0.0346)		(0.0354)	(0.0393)
	0.8 normal 150	-0.0132	-0.0094	-	-0.0204	-0.0179	-
				(0.0482)	(0.0396)		(0.0523)	(0.0439)
			300	-0.0085	-0.0076	-	-0.0114	-0.0103	-
				(0.0312)	(0.0272)		(0.0380)	(0.0333)
			600	-0.0038	-0.0031	-	-0.0040	-0.0037	-
				(0.0217)	(0.0186)		(0.0253)	(0.0216)
	0.8	t(1)	150	-0.0018	-0.0230	-	-0.0025	-0.0192	-
				(0.0243)	(0.0502)		(0.0150)	(0.0465)
			300	-0.0002	-0.0113	-	-0.0006	-0.0113	-
				(0.0061)	(0.0375)		(0.0059)	(0.0308)
			600	0.0001	-0.0077	-	-0.0001	-0.0075	-
				(0.0028)	(0.0263)		(0.0029)	(0.0219)
	0.8	t(3)	150	-0.0138	-0.0151	-	-0.0140	-0.0181	-
				(0.0453)	(0.0455)		(0.0475)	(0.0521)
			300	-0.0054	-0.0073	-	-0.0067	-0.0099	-
				(0.0293)	(0.0338)		(0.0321)	(0.0362)
			600	-0.0038	-0.0052	-	-0.0041	-0.0052	-
				(0.0196)	(0.0227)		(0.0220)	(0.0251)
	(1) In Table 1 for ρ, the biases are small and the standard deviations display a decrease
	pattern as sample size increases. When i ∼ t(1), ρ of the RSSoFR always performs
					20

Table 2 :

 2 The empirical average IMSE and its standard deviations (in brackets) of β(t).

					γ = 1.1			γ = 2
	ρ	i	n	RSSoFR	SSoFR	SoFR	RSSoFR	SSoFR	SoFR
	0	normal 150	0.1386	0.1267	0.1377	0.1727	0.1695	0.1723
				(0.0567)	(0.0531)	(0.0549)	(0.0776)	(0.0755)	(0.0773)
			300	0.0788	0.0721	0.0786	0.1418	0.1403	0.1417
				(0.0346)	(0.0316)	(0.0342)	(0.0455)	(0.0451)	(0.0454)
			600	0.0419	0.0383	0.0418	0.1279	0.1273	0.1279
				(0.0168)	(0.0153)	(0.0168)	(0.0246)	(0.0243)	(0.0247)
	0	t(1)	150	0.2676	30034	0.2632	0.2080	431	0.2084
				(0.1497)	(565196)	(0.1431)	(0.0932)	(6064)	(0.0934)
			300	0.1524	12790	0.1513	0.1569	1356	0.1568
				(0.0752)	(192370)	(0.0745)	(0.0523)	(17250)	(0.0521)
			600	0.0755	4633	0.0754	0.1350	43476	0.1349
				(0.0336)	(71726)	(0.0337)	(0.0267)	(593702)	(0.0266)
	0	t(3)	150	0.1841	0.3278	0.1832	0.1923	0.2204	0.1917
				(0.0885)	(0.2543)	(0.0875)	(0.0859)	(0.1064)	(0.0859)
			300	0.0998	0.1845	0.0996	0.1499	0.1648	0.1500
				(0.0456)	(0.1220)	(0.0452)	(0.0500)	(0.0604)	(0.0501)
			600	0.0535	0.0970	0.0534	0.1290	0.1369	0.1290
				(0.0218)	(0.0700)	(0.0218)	(0.0250)	(0.0327)	(0.0250)
	0.5 normal 150	0.1366	0.1231	0.1984	0.1742	0.1719	0.1946
				(0.0590)	(0.0520)	(0.0904)	(0.0787)	(0.0786)	(0.0882)
			300	0.0761	0.0698	0.1184	0.1410	0.1394	0.1540
				(0.0321)	(0.0281)	(0.0544)	(0.0465)	(0.0461)	(0.0501)
			600	0.0430	0.0395	0.0687	0.1271	0.1265	0.1385
				(0.0171)	(0.0158)	(0.0268)	(0.0255)	(0.0255)	(0.0288)
	0.5	t(1)	150	0.2822	2827	1.1746	0.2104	3744	0.4115
				(0.1471)	(25820)	(0.9887)	(0.0960)	(49109)	(0.3345)
			300	0.1543	3862	0.6101	0.1596	22641	0.2519
				(0.0820)	(49921)	(0.3746)	(0.0525)	(488421)	(0.1293)
			600	0.0740	9897	0.3027	0.1345	134	0.1836
				(0.0335)	(192399)	(0.1699)	(0.0261)	(1001)	(0.0591)
	0.5	t(3)	150	0.1940	0.3148	0.3114	0.1891	0.2216	0.2164
				(0.0929)	(0.1809)	(0.1594)	(0.0861)	(0.1347)	(0.0981)
			300	0.1042	0.1808	0.1722	0.1483	0.1630	0.1676
				(0.0470)	(0.1095)	(0.0826)	(0.0472)	(0.0576)	(0.0560)
			600	0.0531	0.0926	0.0908	0.1289	0.1376	0.1420
				(0.0224)	(0.0591)	(0.0388)	(0.0241)	(0.0329)	(0.0283)
	0.8 normal 150	0.1432	0.1298	0.6378	0.1759	0.1734	0.3524
				(0.0636)	(0.0581)	(0.3425)	(0.0803)	(0.0799)	(0.1792)
			300	0.0763	0.0697	0.4035	0.1425	0.1407	0.2884
				(0.0296)	(0.0264)	(0.1882)	(0.0455)	(0.0450)	(0.1148)
			600	0.0416	0.0377	0.2581	0.1288	0.1282	0.2494
				(0.0168)	(0.0149)	(0.0978)	(0.0248)	(0.0246)	(0.0652)
	0.8	t(1)	150	0.2877	1071	10.7320	0.2029	737	3.0490
				(0.2352)	(7593)	(24.6818)	(0.0916)	(10172)	(9.5959)
			300	0.1518	2446	4.4654	0.1540	932	1.0579
				(0.0730)	(19969)	(4.5632)	(0.0508)	(12959)	(1.0765)
			600	0.0782	7359	1.9833	0.1343	628	0.6228
				(0.0394)	(101428)	(1.3628)	(0.0267)	(7479)	(0.5507)
	0.8	t(3)	150	0.1819	0.3356	0.8924	0.1877	0.2223	0.4125
				(0.0920)	(0.5530)	(0.5139)	(0.0833)	(0.1090)	(0.2617)
			300	0.0989	0.1779	0.5534	0.1470	0.1623	0.3185
				(0.0472)	(0.1357)	(0.2888)	(0.0505)	(0.0561)	(0.1527)
			600	0.0524	0.0899	0.3368	0.1294	0.1362	0.2671
				(0.0216)	(0.0580)	(0.1328)	(0.0256)	(0.0298)	(0.0917)

Table 3 :

 3 The estimates and its standard deviations (in brackets) of ν.

					γ = 1.1			γ = 2
	ρ	i	n	RSSoFR SSoFR	SoFR	RSSoFR SSoFR	SoFR
	0	normal 150	4.0080	-	4.0080	4.0260	-	4.0420
				(0.0892)		(0.1549)	(0.1593)		(0.2105)
			300	4.0040	-	4.0040	4.0000	-	4.0000
				(0.0632)		(0.0632)	(0.0000)		(0.0000)
			600	4.0000	-	4.0000	4.0000	-	4.0000
				(0.0000)		(0.0000)	(0.0000)		(0.0000)
	0	t(1)	150	1.7960	-	1.8260	1.8560	-	1.8820
				(0.4635)		(0.4474)	(0.4557)		(0.4477)
			300	1.8720	-	1.8880	1.8880	-	1.8920
				(0.3462)		(0.3281)	(0.3401)		(0.3355)
			600	1.9560	-	1.9620	1.9620	-	1.9660
				(0.2053)		(0.1914)	(0.1914)		(0.1814)
	0	t(3)	150	3.2300	-	3.2360	3.1660	-	3.2040
				(0.4213)		(0.4250)	(0.3725)		(0.4034)
			300	3.1380	-	3.1440	3.1260	-	3.1580
				(0.3452)		(0.3514)	(0.3322)		(0.3651)
			600	3.0880	-	3.1080	3.0640	-	3.0960
				(0.2836)		(0.3107)	(0.2450)		(0.2949)
	0.5 normal 150	4.0060	-	4.0320	4.0180	-	4.0500
				(0.1342)		(0.1872)	(0.1331)		(0.2601)
			300	4.0020	-	4.0000	4.0000	-	4.0000
				(0.0447)		(0.0000)	(0.0000)		(0.0000)
			600	4.0000	-	4.0000	4.0000	-	4.0000
				(0.0000)		(0.0000)	(0.0000)		(0.0000)
	0.5	t(1)	150	1.7800	-	1.9580	1.8460	-	1.9640
				(0.4817)		(0.6520)	(0.4457)		(0.6321)
			300	1.9100	-	1.9340	1.8840	-	1.8900
				(0.2934)		(0.5120)	(0.3388)		(0.5200)
			600	1.9540	-	1.8900	1.9700	-	1.8600
				(0.2097)		(0.3978)	(0.1708)		(0.3959)
	0.5	t(3)	150	3.2520	-	3.4940	3.1960	-	3.4460
				(0.4346)		(0.5005)	(0.3974)		(0.4976)
			300	3.1360	-	3.5140	3.0900	-	3.3760
				(0.3431)		(0.5003)	(0.2865)		(0.4849)
			600	3.0680	-	3.4600	3.0700	-	3.3920
				(0.2520)		(0.4989)	(0.2554)		(0.4887)
	0.8 normal 150	4.0020	-	4.0540	4.0120	-	4.0520
				(0.0775)		(0.2349)	(0.1260)		(0.2478)
			300	4.0000	-	4.0060	4.0020	-	4.0120
				(0.0000)		(0.0773)	(0.0447)		(0.1090)
			600	4.0000	-	4.0000	4.0000	-	4.0000
				(0.0000)		(0.0000)	(0.0000)		(0.0000)
	0.8	t(1)	150	1.7960	-	2.2380	1.8460	-	2.2240
				(0.4458)		(0.7891)	(0.4457)		(0.7791)
			300	1.8880	-	2.2120	1.9060	-	2.1300
				(0.3157)		(0.6928)	(0.2989)		(0.7114)
			600	1.9540	-	2.0380	1.9660	-	2.0280
				(0.2097)		(0.6398)	(0.1814)		(0.6261)
	0.8	t(3)	150	3.2240	-	3.8120	3.1760	-	3.7780
				(0.4173)		(0.4012)	(0.3812)		(0.4160)
			300	3.1140	-	3.8760	3.0800	-	3.8520
				(0.3181)		(0.3299)	(0.2716)		(0.3555)
			600	3.0820	-	3.9580	3.0580	-	3.8720
				(0.2746)		(0.2008)	(0.2340)		(0.3344)

Table 4 :

 4 The values of the Moran's I statistic under different values of k and d 0

  We choose the W under which the Moran's I statistic takes the greatest value. Set k = {1, 2, 3, 4, 5, 6, 7, 8}. And for each k, we search a d 0 where the value of the Moran's I statistic is the greatest. Table4displays the results for different values of k. It can be seen W with k = 2 and d 0 = 626.0 is the best option. We also found this W performs better than the spatial weight matrix formed by the negative exponential of distance

	km)	613.5	626.0	626.0	626.0	613.5	625.4	625.4	613.5
	Moran's I 0.746	0.801	0.732	0.753	0.754	0.747	0.741	0.743

(2.2).

Table 5 :

 5 The fitting and prediction results of the RSSoFR and the SoFR.

	models	ρ	Moran's I statistic(residuals) MSE(fitted error) MSE(prediction error)	ν
	SoFR	-	0.487	0.128	0.114	-
	RSSoFR 0.675	-0.14	0.088	0.076	15

This research was financially supported by the National Natural Science Foundation of China under grant nos. 71420107025 and 11701023.

similar conclusions as those from results with the FPC basis in Table 1-3. For sake of space, we do not report the corresponding results and will illustrate the FPLS basis in the real data analysis. To sum up, our proposed model performs better than the SoFR and is more robust than the SSoFR, which provides a competitive alternative for the existing methods in FDA.

REAL DATA ANALYSIS

In this section, we revisit the weather dataset presented in Section 1 to assess the application of the RSSoFR. Specifically, we add records that correspond to the weather data for 2008 derived from the China Meteorological Yearbook. In this dataset, records from 2005 to 2007 serve as training set while those from 2008 serve as test set. Note that each of the involved 34 Chinese cities covers a large area with the average land size around 16147 km 2 , and each observation represents the global picture of the monthly mean temperature and monthly total precipitation across the whole city instead of the particular point, which integrates records from tens or even hundreds of meteorological stations in that city area rather than only one weather station (that is commonly treated as point-referenced data). From this perspective, here the dataset can be viewed as areal data instead of point-reference data. We compare the RSSoFR with the SoFR by their fitting and prediction results.

During preprocessing, we smooth the discrete temperature records x i (t j ), i = 1, 2, . . . , 34, t j = 1 12 , 1 6 , . . . , 11 12 , 1 (each j represents a month) using the Epanechnikov Kernel to get x i (t). Figure 3 (left) presents the functional predictors.

Let the response y i be the logarithm of the mean annual total precipitation for the ith city. We build the RSSoFR model as

term has thick tails. An examination of a real dataset illustrates the superiority of the RSSoFR over the SoFR. In short, the robust spatial autoregressive scalar-on-function regression model with t-distribution presented in this paper constitutes a practical statistical tool for modelling the spatial dependent data with functional covariates and scalar response that complements the widely popular spatial autoregressive scalar-onfunction regression with normality assumption and may help to further understanding in many fields of applied research.

We discussed the RSSoFR under the assumption that only one functional predictor is involved. The RSSoFR with multiple functional predictors also deserve attention.

Based on the estimation method introduced in Section 3, the methods presented here can be easily generalized to the RSSoFR with multiple functional predictors.

One referee raised concerns about the model uncertainty. Honestly, apart from the estimation method, the problem of statistical inference is also of great interest as it provides an overall assessment of the association of the functional covariates with the responses (β(t), ρ), as well as the whole model; however, it remains challenging due to the infinite dimensionality of the functional covariates. To overcome this issue, a natural strategy is to reduce the dimension. With representing the functional covariates and the coefficient function by linear combinations of a set of basis functions, the testing problem for β(t) reduces to the hypothesis testing under a classical linear model.

Along with this line, there is a plethora of literature that develops statistical methods.

For example, [START_REF] Su | Hypothesis testing in functional linear models[END_REF] proposed a Wald-type test with varying thresholds in selecting the number of PCs for the functional covariates; [START_REF] Garca-Portugus | A goodness-of-fit test for the functional linear model with scalar response[END_REF] introduced the projected Cramér-von Mises (PcVM) test-a testing method which is derived by using random projection, and whose null distribution is approximated by bootstrap. We refer the reader to [START_REF] Tekbudak | A comparison of testing methods in scalar-on-function regression[END_REF] for an extensive comparison of testing methods in scalar-on-function regression.

However, in the current paper we adopt the EM algorithm to implement the estimation procedure, which cannot obtain the corresponding standard errors of the proposed estimators. A common practice is to compute the standard errors of the proposed estima-