Estimation for spatial semi-functional partial linear regression model with missing response at random

1 Introduction

Nowadays, statistical techniques are essential for the analysis of massive volumes of spatial data, particularly in fields such as environmental sciences, climate, geography, econometrics, medicine, biology, and other applied areas. For example, such data are commonly encountered in environmental and socioeconomic studies, where the goal is to assess the effect of environmental or other ecological parameters on the variability of biomass and the spatial distribution of species or groups of species in marine fauna [1]. These datasets are often sparse, high-dimensional, and include highly correlated responses. This is also the case in oil exploration, where the aim is to predict the physical parameters of oil layers while considering other available parameters from oil fields [2]. Additionally, research on hyperspatial image processing, as developed by [3], involves similar challenges. These data are collected continuously at different stations, and the goal is to understand the relationships between these multiple correlated responses. Consequently, high-dimensional data become functional data, which are then processed and analyzed using functional data analysis (FDA) methods. Such data are referred to as functional data with spatial correlation. Functional geostatistics is a field of application for this type of data. The theoretical and practical developments in this new branch of statistics have led to significant advancements in areas with geographic dependence (for recent reviews on the subject, see [4,5]).

FDA has experienced rapid growth in recent years, evolving from exploratory and descriptive data analysis to more sophisticated linear models and estimation techniques. This dynamic progression has contributed to theoretical and methodological improvements, as well as the expansion of applications across various fields (see, for example, [6–8] and recent bibliographic discussions such as [9,10]). One major area of research in this field is functional regression, which examines the influence of a functional random variable on a scalar variable. Among the key questions in this domain are semi-parametric functional regression models, which combine the flexibility of parametric regression models with the advantages of nonparametric approaches, without the sensitivity to dimensionality issues. For a comprehensive discussion on semi-parametric modeling in functional regression, readers are referred to the bibliographical surveys by [11].

Specifically, one of the most significant semi-parametric functional models is the partially linear regression model introduced by [12], known as the SFPLR model. This model has been widely studied and applied in various contexts. This model is expressed as

where Y is the scalar response variable, X = ( X 1 , X 2 , … , X p ) is a p -dimensional vector of explanatory variables, β is an unknown p -dimensional parameter vector, Z is a functional explanatory variable, m ( ⋅ ) is an unknown smooth functional operator, and ε are identically distributed random errors satisfying E ( ε ) = 0 and unknown variance σ 2 ( ε ) < ∞ . We typically assume an additional condition of independence between the error variable ε and the random vector ( X , Z ) . The authors apply the classical kernel method (Nadaraya-Watson-type weighting method) to establish the asymptotic normality of β and the convergence rate of m . This model was later extended to dependent data by [13]. Furthermore, Aneiros-Pérez and Vieu [14] proposed a bootstrap procedure to approximate the distribution of these estimators, while Lian [15] generalized this model to the case where the linear component is also functional. More recently, Aneiros Pérez et al. [16] examined different bootstrapping procedures for this model under various dependency structures. Additionally, a procedure for testing linearity in partially linear functional models was proposed by Zhao and Zhang [17]. Other methods for estimating the parameters of the SFPLR model include the local linear approach by Feng and Xue [18], robust procedures considered by Boente and Vahnovan [19], the k -nearest neighbors ( k -NN) method used by Ling et al. [20], and Bayesian approaches proposed by Shang [21]. Moreover, Kedir et al. [22] introduced the local linear- k NN approach. For more recent advancements, we refer readers to the bibliographic reviews by previous studies [10,23].

Furthermore, only a few research works have addressed estimation in the SFPLR model for spatially dependent observations. We cite the work of Li and Ying [24] on the autoregressive semi-functional linear regression (SAR) model, which proposed an estimator based on the quasi-maximum-likelihood method and local linear estimation. Additionally, Benallou et al. [25] derived the asymptotic normality of the parametric component, as well as the convergence in probability with the rate of the nonparametric component. The complete analysis of the data is the subject of all the works listed. However, in many applications, this is unfortunately not the case, and we encounter missing data. Common causes of missing responses include faulty equipment, sample contamination, manufacturing defects, clinical study dropouts, climatic conditions, and inaccurate data entry. On the subject of missing data, imputation methods and the uncertainty associated with such imputations have been extensively studied in the multivariate case (see, for example, [26,27]). However, little research has been conducted on nonparametric functional regression with missing data. The first significant work was carried out by Ferraty et al. [28], who estimated the mean of a scalar response using an independent and identically distributed (i.i.d.) functional sample. In this study, the functional regressor is fully observed, but some responses are missing at random (MAR). They generalize the result obtained by Cheng [29] in the multivariate case. Meanwhile, Ling et al. [30] focused on estimating the regression function and established its asymptotic properties under the conditions of stationarity and ergodicity with MAR responses. On the other hand, Rachdi et al. [31] used the k -NN method in combination with the local linear method to estimate the regression function with a small number of randomly missing responses (MAR).

For semi-parametric partially linear multivariate models with randomly missing responses for i.i.d. data, we can cite previous studies [32,33]. The first results obtained for SFPLR models with MAR responses were established by Ling et al. [34], which generalizes those obtained by Wang et al. [32]. Later, Zou et al. [35] considered the partial functional linear model, a special case of the SFPLR model, where the covariate in the parametric component is measured with additive error and the responses are MAR. They used the least-squares method and empirical likelihood methods to construct the model’s estimators.

However, the handling of missing data in spatial data remains underdeveloped. We cite the work of [36] in the multivariate case, which estimates the regression function using the kernel method. The performance of the estimator obtained is compared with the k -NN method. Puranik et al. [37] explored multiple spatial regression models with missing data, using the regression-based imputation method to predict missing values and adding a normally distributed residual term to each predicted value. This approach restores the loss of variability and biases associated with regression imputation. For functional data, Alshahrani et al. [38] proposed studying the kernel estimation of the regression function when the data are spatially functional and MAR. They obtain the asymptotic properties of this estimator, particularly the probability convergence (with rates) and the asymptotic normality of the estimator under certain weak conditions.

In this study, we propose to investigate SFPLR models for spatial data with missing responses. This article is organized as follows: in Section 2, the semi-functional partial linear model for spatial data is presented with these estimators, while the notations and hypotheses used in our work are postponed to Section 3. Section 4 is devoted to the asymptotic results. In Sections 5 and 6, a computer study on simulated and real data is carried out to demonstrate that the proposed estimators allow a nice improvement compared to the usual global approach. Finally, appendix gives some useful lemmas and the detailed proofs of our main results.

2 Model and its estimators

Let Z N be the integer lattice points in the N -dimensional Euclidean space, and let Λ i = { ( Y i , X i ( 1 ) , X i ( 2 ) , … , X i ( p ) , Z i ) T , i = ( i 1 , … , i N ) ∈ Z N } be a measurable strictly stationary spatial process defined over a probability space ( Ω , A , P ) and identically distributed as ( Y , X 1 , X 2 , … , X p , Z ) T = ( Y , X , Z ) T ∈ R × R p × ℱ , where ℱ represents a functional semi-metric space equipped with a semi-metric d . The term site is used to designate a point i . As basic assumption in the non-parametric literature, we assume that the process can be observed in the rectangular domain ℐ n = { i = ( i 1 , … , i N ) ∈ Z N , 1 ≤ i k ≤ n k , k = 1 … , N } , with a sample size of n ^ = n 1 × … × n N , where n = ( n 1 , … , n N ) . Suppose, moreover that, for l = 1 … , N , n l approaches infinity at the same rate: for any j ≠ k ∈ { 1 , … , N } , we have C 1 < ∣ n j ⁄ n k ∣ < C 2 for some 0 < C 1 < C 2 < ∞ , and we write that n → ∞ if min k = 1 … , N ( n k ) → ∞ . This type of observation domain design is known as an asymptotically increasing domain, which allows this domain to expand, while keeping a minimum distance between observation sites. These data whose spatial locations are regular lattices in Z N are the analog of time series observed at equally spaced instances of time. Such data sets are frequently used in practice because the process under study is essentially discrete. To illustrate this last point, it often happens that in order to collect data, by orbiting satellites, on an agricultural surface, which is exploited in certain proportions (for example) by wheat, corn, soybeans, etc. and which must be estimated. These different crops have their reflectance properties, which, together with noise, are detected remotely. Then, the procedure is to divide the surface of this Earth into small rectangles called pixels. Thus, the data are received as a regular lattice in Z 2 and are identified with the centers of their respective pixels.

The main goal is to predict the spatial process ( Y i , i ∈ Z N ) in some unobserved locations, particularly at an unobserved site i 0 ∈ ℐ n under the information that can be drawn on ( X , Z ) i 0 and observations ( Λ i ) i ∈ ϑ , where ϑ is the observed spatial set of finite cardinality tending to ∞ as n → + ∞ and contained in ℐ n , with i 0 ∉ ϑ .

As mentioned in Section 1, to meet this objective, we consider, in this work, the case where the response variable Y depends on a p -dimensional random vector X in a linear way but is related to functional variable Z in an indeterminate form through an SFPLR model:

where X i = ( X i ( 1 ) , X i ( 2 ) , … , X i ( p ) ) , β = ( β 1 , … , β p ) T , m ( ⋅ ) and ε are defined as before such that

We denote by B ( Z , h ) = { Z ′ ∈ ℱ such that d ( Z , Z ′ ) ≤ h } the topological closed ball.

The kernel estimators of β n and m n (see [25]) are defined by

and

where

with

where K denotes a real-valued kernel function and b n a decreasing sequence of bandwidths, which tends to zero as n tends to infinity.

In this article, we develop an estimation approach in which the values of the independent variables ( X and Z ) are all observed; however, some observations of the response variable ( Y values) are missing. We recall that we say that Y i is missing if Y i does not contain all of the required elements.

For this, we consider a real random variable denoted δ such that δ i = 1 if the value Y i is known, and δ i = 0 , otherwise. Thus, the study will be carried out on an incomplete sample of size n :

It is assumed that the random missing data mechanism satisfies the saving condition:

This conditional probability p ( x , z ) is generally unknown.

First, we note that

From our assumption, by conditioning on Z i = z , it follows that

from where we have

In the following, we define m X ( z ) = E ( δ X ∣ Z = z ) ⁄ E ( δ ∣ Z = z ) , m Y ( z ) = E ( δ Y ∣ Z = z ) ⁄ E ( δ ∣ Z = z ) , and we write A ⊗ 2 = A A T .

Then, by equations (5) and (6), we can write

Thus, if the functions m X ( Z ) and m Y ( Z ) are known, the least-squares estimator of β is given by

However, m X ( Z i ) and m Y ( Z i ) are generally unknown, and they must be estimated for using 8. Assuming that m X and m y are smooth functions of Z i , these can be estimated by using the nonparametric Nadaraya-Watson kernel estimator noted by m ˆ X ( Z ) and m ˆ Y ( Z ) and expressed by

where

with K being a real-valued kernel function and h n a decreasing sequence of bandwidths, which tends to zero as n tends to infinity.

Hence, an estimator of β n is given by

where

and

We then insert β ˆ n in equation (7) to obtain

3 Notations and hypotheses

Our main objective is to obtain the asymptotic normality of our estimator under the condition that the process Λ i is strictly stationary, which satisfies the following α -mixing condition [39]: there exists a real function φ ( t ) , which tends to 0 when t goes to ∞ , such that for finite cardinal subsets (Card ( ⋅ ) < ∞ ) E , E ′ ⊂ Z N :

where d ′ ( E , E ′ ) denotes the Euclidean distance in Z N , B ( S ) = B ( Z i , i ∈ S ) , for S ⊂ Z N , the σ -fields generated by a random variable Z i , and ψ : Z 2 ⟶ R + is a symmetric positive function nondecreasing in each variable such that

Moreover, as is often the case in spatial regression, we assume other conditions on the function φ :

Noting that both conditions (12) and (13) are satisfied by many spatial models (see, for example, the spatial linear process [40,41] for more details on the required conditions). It should be pointed out that if ψ = 1 , then the process Λ i is called a strong mixing process.

In the following, and in order to simplify the reading of the manuscript, we introduce the following notations.

Let ν ij = sup i ≠ j P [ ( Z i , Z j ) ∈ ( B ( z , h ) × B ( z , h ) ) ] the probabilistic joint distribution of Z i and Z j , and let ϕ z ( h n ) = P ( Z ∈ B ( z , h n ) ) called the small ball probability.

In addition, we note θ i ( s ) ( z ) = X i ( s ) − E [ δ i X i ( s ) ∣ Z i = z ] , for s = 1 , … , p , and θ i = ( θ i ( 1 ) , … , θ i ( p ) ) T . Recall that the expressions of our estimators in (9) and (10) contain estimators of θ i . Now, in order to state the asymptotic properties of our estimators, we introduce the following assumptions.

(H1) We suppose that, ∀ z ∈ ℰ and ∀ i ≠ j ∈ Z N , there exist C > 0 such that

for some 1 < a < γ N − 1 , sup i ≠ j ν ij ≤ C ( ϕ z ( h ) ) a + 1 a .

(H2) The kernel function K : R ⟶ R + is assumed to be differentiable with support in the interval [ 0 , 1 ] such that there exist two constants C 3 and C 4 with

(H3) We assume that any function f ∈ { m , θ i ( 1 ) ( z ) , … , θ i ( p ) ( z ) } is smooth, i.e., there exists c > 0 such that for κ > 0 , we have

(H4) We suppose that, ∀ z ∈ ℰ : ϕ z ( h n ) > 0 , there exists differentiable nonnegative functions τ and g such that

(H5) For all t ∈ [ 0 , 1 ] , we have

(H6) Let R i = θ i ε i and Σ = E [ p ( X 1 , Z 1 ) ( θ 1 ( θ 1 ) T ) ] , where 1 is the site spatial ( 1 , … , 1 ) and we denote by 0 the site spatial ( 0 , … , 0 ) .

1. The matrix B = ∑ i ∈ ℐ n E [ p ( X i , Z i ) ( R 0 ( R i ) T ) ] is assumed to be positive definite.

2. We assume that Σ is an invertible matrix.

(H7) We suppose that

1. E ∣ ε 1 ∣ ρ + E ∣ θ 1 ( 1 ) ∣ ρ + … + E ∣ θ 1 ( p ) ∣ ρ < ∞ , for some ρ ≥ 3 .

2. For all i ≠ j , E [ Y i Y j ∣ ( Z i , Z j ) ] < ∞ .

3. For all i ≠ j , max 1 ≤ s ≤ p E [ X i ( s ) X j ( s ) ∣ ( Z i , Z j ) ] < ∞ .

(H8) Let V 2 ( z ) = v a r ( Y i − X i T β ∣ Z i = z ) and V s ( z ′ ) = E [ ∣ Y i − X i T β − m ( z ) ∣ s ∣ Z i = z ′ ] , with s > 2 .

1. We suppose that the functions V 2 ( ⋅ ) and V s ( ⋅ ) are continuous functions near z , i.e., since h → 0 , we have

   sup { z ′ : d ( z , z ′ ) ≤ h } ∣ V k ( z ′ ) − V k ( z ) ∣ = o ( 1 ) , k ≥ 2 .

2. Let V ( z ′ , z ″ , z ) = E [ ( Y i − X i T β − m ( z ) ) ( Y j − X j T β − m ( z ) ) ∣ Z i = z ′ , Z i = z ″ ] , for i ≠ j .

(H9) We also suppose that p 1 ( ⋅ ) = P ( δ = 1 ∣ Z = z ) is continuous function near z , i.e.,

(H10) There exists 1 > ϑ > N γ such that ϕ z ( h ) ≥ n ^ ϑ − 1 2 N + 1 .

Comments on the assumptions:

The asymptotic results obtained subsequently are based on Theorem 4.1 in [38]. It is, therefore, normal to impose the same hypotheses as in this theorem. More specifically, hypothesis ( H 1 ) measures the local dependence between observations, and such condition is necessary to reach the same rate of convergence as in the i . i . d . case. The assumption ( H 2 ) is also standard in this context of functional statistics. The conditions in these assumptions are satisfied by several kernels, which are traditionally used in practice such as Epanechnikov kernel, Parzen kernels, triangular kernel, and others. Hypothesis ( H 2 ) and hypothesis ( H 3 ) are used to control the bias of the estimator. The assumptions ( H 4 ) is known as the “concentration property” in the asymptotic theory of nonparametric functional statistics. It is a condition closely linked with the semi-metric d and such function ϕ z ( ⋅ ) can be explicated for several continuous processes [8]. The information on the variability of the small-ball probability is given by hypothesis ( H 5 ) . This information is used to control the bias of the nonparametric estimators. Moreover, hypothesis ( H 5 ) is necessary to apply the central limit theorem to have asymptotic normality. It is well documented that this condition is not very restrictive and the function ι z ( t ) can be given explicitly for many classical stochastic processes [42]. Assumptions ( H 8 ) and ( H 10 ) present the continuous local conditions needed to establish the main results. In fact, the conditional expectation requires that V s and V are continuous. Assumption ( H 9 ) in FDA MAR models is typical (see, for example, [30,31]). The additional hypothesis ( H 6 ) is usual in the context of the SFPLR.

4 Theoretical results

We are now in a position to give our asymptotic results. The first one gives the asymptotic distribution of the estimator for the parametric component of the model ( β ^ n ).

The following results give the probability convergence and the asymptotic normality of the estimator of the non-parametric part.

We note that these results obtained extend those that are established in the case of complete data [25].

5 Computational study

In this section, we are interested in the behavior of the estimators proposed on samples of finite size, with particular attention to the influence of spatial correlation and the effect of MAR on the efficiency of the estimators. To do that, we conducted simulations based on observations denoted as ( X i , Z i , Y i , δ i ) with i = ( i 1 , i 2 ) , 1 ≤ i 1 ≤ n 1 , 1 ≤ i 2 ≤ n 2 , and we generate the SFPLR model with MAR as follows:

where β = ( 1 , 2 ) T , X i j ∽ Exponential ( 0.5 ) , j = 1 , 2 , and take the nonparametric operator m as follows:

The curves Z i ( t ) are defined as follows: Z i ( t ) = B i t sin ( t − A i ) + ν i ( t ) , for t ∈ [ 0 , 1 ] .

For simulating the curves Z i , we take ν i ( t ) ∽ N ( 0 , 0.2 ) , A = D * sin G 2 + 0.5 , with G = G R F ( 0 , 5 , 3 ) , B = G R F ( 2.5 , 5 , 3 ) , and ε = G R F ( 0 , 0.1 , 5 ) , where the function G R F ( μ , σ 2 , s ) denotes a stationary Gaussian random field with mean μ and covariance function defined by C ( l ) = σ 2 exp − ‖ l ‖ s 2 , l ∈ R 2 and s > 0 . And the function D is defined by

The curves, following the values of a , are displayed in Figure 1. Note that we consider the same curves as in [43], where the function D is used to control the spatial mixing condition. Therefore, these observations are a mixture of independent and dependent data points, as illustrated in Figure 2 (for more details, see [43]).

Figure 1

Curves Z i , t ∈ [ 0 , 1 ] for a = 5 , 20, and 50.

Figure 2

Simulations of the random field were generated for different values of a , specifically a = 5 , 20, and 50.

To reduce the level of independence in the data, it is sufficient to decrease the value of a . For our analysis, we have employed a value of σ 2 = 5 and s = 5 . Moreover, similar to that described in [34], we adopted the following missing data mechanism:

where expit ( u ) = e u ⁄ ( 1 + e u ) , for all u ∈ R , and we will take for α the following values: α = 0.05 , 0.5, and 5.

Let us remember that the degree of dependence between the functional variable Z and the variable δ is controlled by the parameter α , and to check the value of p ( x ) , we compute δ ¯ = 1 − 1 n 1 × n 2 ∑ i 1 = 1 n 1 ∑ i 2 = 1 n 2 δ ( i 1 , i 2 ) for different values of α to measure the missing rate.

In order to calculate our estimators, we use the class of semi-metrics based on the principal component analysis (PCA) metric, which is best suited to the treatment of this type of data (discontinuous functional variable). Furthermore, we have chosen the standard kernel function defined as follows: K ( u ) = 3 2 ( 1 − u 2 ) 1 [ 0 , 1 ] ( u ) .

The objective of this computer study is to conduct a comparison between our semi-partially linear model with MAR estimator (SFPLRM), semi-partially linear model estimator r ˜ n in the complete case (SFPLRC) (see [25]), and the nonparametric functional model estimator with MAR (FNPM).

Recall that the FNPM model is defined as

and the estimator of operator regression r [38] is given by

For that purpose, we conducted a random split of our data, denoted as ( X i , Z i , Y i , δ i ) i , into two subsets: a test sample ( X i , Z i , Y i ) i ∈ I ′ (consisting of data without missing values) and training sample ( X i , Z i , Y i , δ i ) i ∈ I that they will be used to select the optimal smoothing parameters h opt . The training sample was used to choose the smoothing parameter ( h opt ). The h o p t is the data-driven bandwidth obtained by a cross-validation procedure:

with m ^ n ( − i ) ( ⋅ ) being the estimator of m ( ⋅ ) based on the leave-one-out method calculated without observation ( X i , Z i ) after estimating β n (see [14,44] for more details).

To evaluate the precision of the estimators of the three models (SFPLRM, SFPLRC, and FNPM), we use the square and mean square errors, denoted mean squared error (MSE):

where # ( I ′ ) represents the size of the testing sample I ′ . The experiment was replicated M = 100 times, which allows us to compute M values for MSE and display their distribution through a boxplot. The results obtained (the prediction values compared with the true values) are presented in Figure 3 for the three models with different values of a . Figure 4 displays the boxplots constructed from MSE obtained for the three models.

Figure 3

Predictions of the three models for a = 5 , 20 , and 50.

Figure 4

Boxplot MSE of the three models for α = 0.05 , 0.5, and 5.

In Figure 4, we observe that the SFPLRM estimator shows better prediction effects than the FNPM in comparison with the SFPLRC estimator.

In Table 1, which presents the MSE for FNPM, SFPLRM, and SFPLRC under various conditions defined by combinations of n 1 , n 2 , and α , it becomes evident that the SFPLRM estimator consistently shows superior predictive accuracy when compared to the FNPM estimator across a range of settings. This is notably evident as the MSE values for the SFPLRM estimator consistently remain lower across the majority of combinations of n 1 , n 2 , and α . Furthermore, a notable observation is that higher values of α (e.g., α = 5 ) tend to correspond with lower MSE values, indicating enhanced predictive performance for the SFPLRM estimator. It is also worth mentioning that as the number of samples ( n 1 × n 2 ) increases, the results suggest that the MAR mechanism has little to no discernible effect on the prediction of MSE. This observation implies that the SFPLRM estimator maintains consistent performance regardless of the presence of missing data, mainly when dealing with larger sample sizes.

Now, we aim to examine the bias in the estimation of β and its associated square error (ASE), which is defined as

The results are presented in Table 2.