A Variance–Covariance Framework for Multiple Nonparametric Regression with Missing Data

1. Introduction

Nonparametric regression has become an indispensable tool for analyzing complex data structures in which the relationship between a response variable and multiple predictors cannot be adequately described by linear or low-dimensional parametric models. Foundational monographs and surveys such as [1,2,3] establish kernel smoothing and local polynomial regression as flexible alternatives to parametric modeling, with wide-ranging applications in economics, biostatistics, environmental sciences, and engineering. Classical developments in this area, however, typically assume that both covariates and responses are fully observed.

In many modern applications, incomplete observations arise naturally due to survey nonresponse, censoring in survival and reliability studies, attrition in longitudinal data, or technical limitations in data collection. The statistical analysis of missing data has therefore developed into a major research area in its own rigfan2018, with comprehensive treatments provided by [4,5]. When missingness is present in nonparametric regression, the interaction between smoothing methods and incomplete data mechanisms creates challenges that are substantially more intricate than in parametric models.

Missing data introduce two fundamental difficulties in nonparametric regression. First, the effective sample size available for local smoothing is reduced, which directly inflates the variance of kernel and local polynomial estimators [1]. Second, and more critically, missingness alters the joint distribution of the covariates and the response, thereby distorting the local structure upon which nonparametric estimators rely. Since smoothing methods are inherently local, even mild distortions of joint moments can lead to substantial bias and instability.

A substantial body of literature has addressed bias correction and consistency under missing data through inverse probability weighting, imputation, and augmented or doubly robust estimating equations; e.g., [6,7,8]. Extensions of these ideas to semiparametric and nonparametric regression settings have been explored by several authors, often focusing on mean regression or partially linear structures. Nevertheless, much of this work treats the covariance structure implicitly, without examining how missingness reshapes the second-order geometry that underlies local smoothing.

The central objective of this paper is to articulate the conceptual and theoretical role of the variance–covariance matrix in multiple nonparametric regression with missing data. We argue that this matrix, although it does not parameterize the regression function itself, acts as a latent geometric object that governs local identifiability, estimator stability, efficiency, and the validity of asymptotic inference. This role becomes especially pronounced under missing at random (MAR) and missing not at random (NMAR) mechanisms, where the observed covariance structure differs systematically from its population counterpart [4].

By emphasizing covariance geometry rather than specific estimators, this paper provides a unifying perspective on why nonparametric regression procedures may fail or become unstable under missing data, and why principled correction methods must account not only for bias but also for distortion in local second-order structure. The discussion that follows aims to bridge the gap between classical nonparametric smoothing theory and modern missing data methodology.

2. Model Framework

Let $(Y,\mathbf{X})$ denote a scalar response variable and a p-dimensional covariate vector, with the multiple nonparametric regression model

$$Y=m\left(\mathbf{X}\right)+\epsilon ,$$

(2.1)

where $m(·):{\mathbb{R}}^p\to \mathbb{R}$ is an unknown smooth regression function capturing the potentially nonlinear relationship between the covariates and the response, and $\epsilon$ is a random error term satisfying

$$\mathbb{E}(\epsilon \mid \mathbf{X})=0\mathrm{and}Var(\epsilon \mid \mathbf{X})={\sigma }^2\left(\mathbf{X}\right)<\infty .$$

The conditional mean function $m\left(\mathbf{X}\right)$ represents the target of inference, while the conditional variance function ${\sigma }^2\left(\mathbf{X}\right)$ allows for heteroscedasticity, which is often encountered in practical applications such as economics, biostatistics, and environmental studies.

Define the joint vector $\mathbf{Z}={(Y,{\mathbf{X}}^{\top })}^{\top }\in {\mathbb{R}}^{p+1}$ and its population variance–covariance matrix

$$\Sigma =Var\left(\mathbf{Z}\right)=\left(\begin{array}{cc}Var\left(Y\right)& Cov(Y,{\mathbf{X}}^{\top })\\ Cov(\mathbf{X},Y)& Var\left(\mathbf{X}\right)\end{array}\right).$$

(2.2)

Here, $Var\left(\mathbf{X}\right)$ is a $p\times p$ positive semidefinite matrix capturing the marginal variation and linear dependence structure among the covariates, while $Cov(Y,{\mathbf{X}}^{\top })$ is a $1\times p$ vector of covariances between the response and each covariate. Although $m\left(\mathbf{X}\right)$ is generally nonlinear and not determined solely by second-order moments, the covariance structure plays a crucial role in local polynomial and kernel-based nonparametric estimation, as it affects identifiability, stability, and efficiency of estimators.

In realistic settings, components of $(Y,\mathbf{X})$ may be incompletely observed. Let $R=(R_Y,R_1,\dots ,R_p)$ denote a vector of missingness indicators, where $R_Y=1$ if Y is observed and 0 otherwise, and $R_j=1$ if $X_j$ is observed and 0 otherwise for $j=1,\dots ,p$. The missingness mechanism is formally defined by the conditional distribution $P(R\mid Y,\mathbf{X}),$ which leads to three classical types of missingness [4]:

(i)

Missing Completely at Random (MCAR):$P(R\mid Y,\mathbf{X})=P\left(R\right)$; missingness is independent of both observed and unobserved data.

(ii)

Missing at Random (MAR):$P(R\mid Y,\mathbf{X})=P(R\mid {\mathbf{X}}_{\mathrm{obs}})$; missingness depends only on observed data.

(iii)

Missing Not at Random (NMAR):$P(R\mid Y,\mathbf{X})$ depends on unobserved values, requiring explicit modeling to ensure identifiability.

The presence of missingness affects both the observed covariance structure and the effective sample size, which in turn influences local nonparametric smoothing, estimator variance, and asymptotic properties. Properly accounting for the missingness mechanism is therefore essential to obtain consistent, efficient, and unbiased estimates of the regression function $m\left(\mathbf{X}\right)$.

3. Variance–Covariance Structure in Multiple Nonparametric Regression

Consider the multiple nonparametric regression model

$$Y=m\left(\mathbf{X}\right)+\epsilon ,\mathbb{E}(\epsilon \mid \mathbf{X})=0,Var(\epsilon \mid \mathbf{X})={\sigma }^2\left(\mathbf{X}\right),$$

(3.1)

where $Y\in \mathbb{R}$ is the response, $\mathbf{X}\in {\mathbb{R}}^p$ is the covariate vector, and $m(·)$ is an unknown smooth regression function. Unlike linear regression, where Gaussian assumptions imply that second-order moments fully characterize the regression function, $m(·)$ in nonparametric settings is generally nonlinear and not determined solely by covariance structure. Nevertheless, the variance–covariance geometry of $\mathbf{X}$ plays a fundamental role in the feasibility, stability, and efficiency of nonparametric estimation.

3.1. Local Polynomial and Kernel Estimation

Kernel and local polynomial methods [1,2] estimate $m\left(\mathbf{x}\right)$ by fitting a local Taylor expansion around a target point $\mathbf{x}$:

$$m\left({\mathbf{X}}_i\right)\approx m\left(\mathbf{x}\right)+{({\mathbf{X}}_i-\mathbf{x})}^{\top }\nabla m\left(\mathbf{x}\right)+\dots ,$$

(3.2)

where $\nabla m\left(\mathbf{x}\right)$ denotes the gradient vector. The corresponding local design matrix is

$$\mathbf{U}\left(\mathbf{x}\right)=\left(\begin{array}{ccc}1& {({\mathbf{X}}_1-\mathbf{x})}^{\top }& \dots \\ ⋮& ⋮& \ddots \\ 1& {({\mathbf{X}}_n-\mathbf{x})}^{\top }& \dots \end{array}\right),$$

(3.3)

and the weighted least squares solution relies on the weighted moment matrix

$$\mathbf{M}\left(\mathbf{x}\right)={\displaystyle \sum _{i=1}^n}{K}_h({\mathbf{X}}_i-\mathbf{x}){\mathbf{U}}_i\left(\mathbf{x}\right){\mathbf{U}}_i{\left(\mathbf{x}\right)}^{\top },$$

(3.4)

where $K_h(·)$ is a kernel function with bandwidth h. The invertibility and conditioning of $\mathbf{M}\left(\mathbf{x}\right)$ depend directly on the local variation of the covariates.

3.2. Local Covariance Matrix

A central object is the local covariance matrix of $\mathbf{X}$ at $\mathbf{x}$:

$${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)=\mathbb{E}\left[(\mathbf{X}-\mathbf{x}){(\mathbf{X}-\mathbf{x})}^{\top }|\mathbf{X}\approx \mathbf{x}\right].$$

(3.5)

This matrix captures the local dispersion and directional dependencies of the covariates. Positive definiteness ensures that sufficient independent variation exists along all directions of $\mathbf{X}$, enabling stable local polynomial fits and reliable kernel smoothing [3]. Near-singularity or ill-conditioning of ${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)$ increases estimator variance and may lead to convergence issues, particularly in regions with sparse observations or strong multicollinearity.

3.3. Implications for Bias–Variance Trade-Off

The local covariance structure of covariates plays a central role in controlling the bias and variance of nonparametric estimators. Consider the local linear estimator of $m\left(\mathbf{x}\right)$:

$$\widehat{m}\left(\mathbf{x}\right)={\mathbf{e}}_1^{\top }{\left(\mathbf{U}{\left(\mathbf{x}\right)}^{\top }\mathbf{W}\left(\mathbf{x}\right)\mathbf{U}\left(\mathbf{x}\right)\right)}^{-1}\mathbf{U}{\left(\mathbf{x}\right)}^{\top }\mathbf{W}\left(\mathbf{x}\right)\mathbf{Y},$$

(3.6)

where $\mathbf{U}\left(\mathbf{x}\right)$ is the local design matrix (including an intercept and linear terms), $\mathbf{W}\left(\mathbf{x}\right)=diag\left\{K_h({\mathbf{X}}_i-\mathbf{x})\right\}$ is the kernel weight matrix, and ${\mathbf{e}}_1$ selects the intercept term corresponding to $m\left(\mathbf{x}\right)$.

Define the weighted local moment matrix:

$$\mathbf{M}\left(\mathbf{x}\right)=\mathbf{U}{\left(\mathbf{x}\right)}^{\top }\mathbf{W}\left(\mathbf{x}\right)\mathbf{U}\left(\mathbf{x}\right)\approx n{f}_{\mathbf{X}}\left(\mathbf{x}\right){\Sigma }_{\mathbf{U}}\left(\mathbf{x}\right),$$

(3.7)

where $f_{\mathbf{X}}\left(\mathbf{x}\right)$ is the marginal density of $\mathbf{X}$ at $\mathbf{x}$, and ${\Sigma }_{\mathbf{U}}\left(\mathbf{x}\right)$ is the second-moment matrix of covariates in the local neighborhood. Its invertibility ensures local identifiability and numerical stability of $\widehat{m}\left(\mathbf{x}\right)$.

The asymptotic bias of $\widehat{m}\left(\mathbf{x}\right)$ for a p-dimensional covariate vector under standard smoothness assumptions is approximately:

$$Bias\left[\widehat{m}\left(\mathbf{x}\right)\right]\approx {\displaystyle \frac{1}{2}}{\displaystyle \sum _{j,k=1}^p}{\displaystyle \frac{{\partial }^{2}m\left(\mathbf{x}\right)}{\partial x_j\partial x_k}}{\mu }_{2,jk}\left(K\right)h^2,$$

(3.8)

where ${\mu }_{2,jk}\left(K\right)=\int u_j{u}_{k}K\left(u\right)du$ are the second moments of the kernel. The local covariance matrix enters directly because ${\mu }_{2,jk}\left(K\right)$ is effectively weighted by the local variation of $\mathbf{X}$ around $\mathbf{x}$; small variance along a direction j reduces the effective second moment, increasing the bias along that direction if the bandwidth h is not appropriately adjusted.

The conditional variance of $\widehat{m}\left(\mathbf{x}\right)$ is

$$Var[\widehat{m}\left(\mathbf{x}\right)\mid \mathbf{X}]\approx {\displaystyle \frac{{\sigma }^2\left(\mathbf{x}\right)}{n{h}^p{f}_{\mathbf{X}}\left(\mathbf{x}\right)}}{\mathbf{e}}_1^{\top }{\Sigma }_{\mathbf{U}}{\left(\mathbf{x}\right)}^{-1}\mathbf{R}\left(K\right){\Sigma }_{\mathbf{U}}{\left(\mathbf{x}\right)}^{-1}{\mathbf{e}}_1,$$

(3.9)

where $\mathbf{R}\left(K\right)=\int K{\left(u\right)}^{2}du$. Ill-conditioned ${\Sigma }_{\mathbf{U}}\left(\mathbf{x}\right)$ inflates variance, confirming that regions with low local variability require larger bandwidths to stabilize estimation. Conversely, well-conditioned local covariance allows smaller bandwidths and finer resolution, reducing bias without substantial variance inflation.

Thus, the bias–variance trade-off is intimately tied to the geometry of ${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)$. Explicitly accounting for its structure enables:

• Adaptive bandwidth selection proportional to local variance.
• Covariance-aware diagnostics to identify regions of potential estimator instability.
• Robust estimation under data irregularities or missingness that distort local dispersion.

Formally, the local covariance matrix acts as a latent geometric object, governing the feasibility, stability, and efficiency of nonparametric regression, providing a rigorous basis for both theoretical and practical implementation.

3.4. Illustrative Numerical Example

To illustrate the practical relevance of the local covariance structure, we consider a small numerical example. Suppose we have $n=6$ observations with $p=2$ covariates and a response variable:

$$\mathbf{X}=\left(\begin{array}{cc}0.1& 0.3\\ 0.2& 0.5\\ 0.4& 0.2\\ 0.6& 0.7\\ 0.8& 0.9\\ 1.0& 0.4\end{array}\right),Y=\left(\begin{array}{c}1.2\\ 1.5\\ 1.1\\ 1.8\\ 2.1\\ 1.7\end{array}\right).$$

For a target point ${\mathbf{x}}_0={(0.5,0.5)}^{\top }$, we compute the local covariance matrix using a simple uniform kernel selecting points within a radius $h=0.3$:

$${\Sigma }_{\mathbf{X}}\left({\mathbf{x}}_0\right)={\displaystyle \frac{1}{\left|N\right({\mathbf{x}}_0\left)\right|}}\sum _{i\in N\left({\mathbf{x}}_0\right)}({\mathbf{X}}_i-{\mathbf{x}}_0){({\mathbf{X}}_i-{\mathbf{x}}_0)}^{\top },$$

where $N\left({\mathbf{x}}_0\right)$ denotes the indices of points within the bandwidth. In this case, the neighbors are observations $i=2,3,4$, giving

$${\Sigma }_{\mathbf{X}}\left({\mathbf{x}}_0\right)=\left(\begin{array}{cc}0.0133& 0.0033\\ 0.0033& 0.0267\end{array}\right).$$

The condition number of this matrix is approximately 2.0, indicating a well-conditioned local design. This suggests that local polynomial or kernel regression around ${\mathbf{x}}_0$ can be performed reliably.

If we consider another target point ${\mathbf{x}}_1={(0.9,0.8)}^{\top }$, the local covariance matrix computed from neighbors $i=5,6$ is

$${\Sigma }_{\mathbf{X}}\left({\mathbf{x}}_1\right)=\left(\begin{array}{cc}0.005& -0.01\\ -0.01& 0.005\end{array}\right),$$

with a condition number exceeding 10, indicating near-singularity. This demonstrates that in regions with sparse observations or highly collinear covariates, the local covariance structure can lead to unstable estimation, highlighting the practical importance of monitoring ${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)$ for robust nonparametric regression.

4. Effects of Missing Data on Covariance Geometry

Missing data can substantially distort the covariance structure of the covariates and response, with implications that are highly dependent on the underlying missingness mechanism [4,5]. Let $\mathbf{Z}={(Y,{\mathbf{X}}^{\top })}^{\top }$ denote the full data vector and R the corresponding missingness indicator. The observed covariance matrix is computed from the available data and may differ from the population covariance matrix $\Sigma =Var\left(\mathbf{Z}\right)$.

4.1. Missing Completely at Random (MCAR)

Under MCAR, the probability of missingness is independent of both observed and unobserved data, i.e., $P(R\mid Y,\mathbf{X})=P\left(R\right)$. In this case, the observed covariance matrix remains an unbiased estimator of the population covariance, although the reduction in sample size increases the variance of the estimators. Local covariance matrices ${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)$ used in kernel or local polynomial smoothing are thus consistent but noisier, potentially requiring larger bandwidths to stabilize the estimates.

4.2. Missing at Random (MAR)

Under MAR, the missingness of a variable depends only on observed data: $P(R\mid Y,\mathbf{X})=P(R\mid {\mathbf{X}}_{\mathrm{obs}})$. Conditioning on observed covariates can introduce systematic distortions in the joint moments of $(Y,\mathbf{X})$, as the observed data no longer represent a simple random sample from the full distribution. This leads to biased or inefficient estimates if the missingness mechanism is ignored. Correction strategies, such as inverse probability weighting (IPW) or augmented estimating equations, are required to recover unbiased estimates of covariance structure and to maintain identifiability of the regression function [6,8].

4.3. Missing Not at Random (NMAR)

Under NMAR, the missingness depends on unobserved values of the data itself, i.e., $P(R\mid Y,\mathbf{X})$ involves unobserved components. In this scenario, the observed covariance structure may fundamentally differ from the population covariance, and standard nonparametric estimators are generally biased without additional modeling assumptions. Estimation requires explicit modeling of the missingness mechanism or the joint distribution of $(Y,\mathbf{X},R)$ to ensure identifiability.

4.4. Implications for Local Covariance Matrices

The impact of missingness propagates directly to the local covariance matrices ${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)$ employed in nonparametric smoothing. In regions where missingness is substantial or systematically biased, the local covariance may become ill-conditioned or nearly singular, undermining the stability of kernel or local polynomial estimators. Consequently, missing data can lead to both increased variance and potential bias in estimated regression surfaces, highlighting the need to account for the covariance distortion induced by incomplete observations during both estimation and bandwidth selection.

5. Kernel and Local Polynomial Estimation

To estimate the regression function $m\left(\mathbf{x}\right)$ at a target point $\mathbf{x}\in {\mathbb{R}}^p$, we consider a local polynomial approach of order $q\ge 0$. Let $K_h(·)=h^{-p}K(·/h)$ denote a multivariate kernel function with bandwidth $h>0$, where $K(·)$ is a symmetric probability density function. The local polynomial estimator is obtained by solving

$$\widehat{\beta }\left(\mathbf{x}\right)=arg\underset{\beta }{min}{\displaystyle \sum _{i=1}^n}{K}_h({\mathbf{X}}_i-\mathbf{x}){\left(Y_i-{\beta }^{\top }{\mathbf{U}}_i\left(\mathbf{x}\right)\right)}^2,$$

(5.1)

where ${\mathbf{U}}_i\left(\mathbf{x}\right)={(1,{({\mathbf{X}}_i-\mathbf{x})}^{\top },\dots ,{({\mathbf{X}}_i-\mathbf{x})}^{\otimes q})}^{\top }$ collects polynomial terms up to degree q centered at $\mathbf{x}$. The estimator of the regression function at $\mathbf{x}$ is then given by

$$\widehat{m}\left(\mathbf{x}\right)={\mathbf{e}}_1^{\top }\widehat{\beta }\left(\mathbf{x}\right),$$

where ${\mathbf{e}}_1$ selects the intercept term.

5.1. Asymptotic Covariance and Local Covariance Structure

The asymptotic covariance of $\widehat{m}\left(\mathbf{x}\right)$ depends on the localized moment matrix

$$\mathbf{M}\left(\mathbf{x}\right)=\mathbb{E}\left[K_h(\mathbf{X}-\mathbf{x})\mathbf{U}\left(\mathbf{x}\right)\mathbf{U}{\left(\mathbf{x}\right)}^{\top }\right],$$

(5.2)

which aggregates the weighted cross-products of the covariate polynomial terms around $\mathbf{x}$. This matrix encodes the local second-order geometry of the covariates and is closely related to the local covariance matrix ${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)$ discussed in Section 4. Its invertibility and condition number directly affect the stability and variance of the local polynomial estimator [2,3].

5.2. Impact of Missing Data

When missing values are present, the observed design matrix is constructed from the available data only, altering both $\mathbf{M}\left(\mathbf{x}\right)$ and the effective sample size. Under MCAR, the matrix remains unbiased but becomes noisier due to reduced observations. Under MAR, the weighting induced by the missingness mechanism modifies the local moments, potentially introducing bias if uncorrected. Under NMAR, the observed local moments may deviate systematically from the population quantities, compromising identifiability and leading to biased or unstable local fits.

Thus, missing data propagate through the local covariance geometry to the estimation step: ill-conditioned $\mathbf{M}\left(\mathbf{x}\right)$ leads to inflated variance or numerical instability, and distorted joint moments may induce systematic bias. Proper handling via weighting, imputation, or augmented estimation is essential to maintain both stability and consistency of the nonparametric estimator in the presence of incomplete data.

6. Simulation Study

To empirically illustrate the effects of missing data on multiple nonparametric regression, we conducted a comprehensive simulation study emphasizing the role of variance–covariance structure. The goal was to examine how different missing data mechanisms affect estimator bias, variance, and stability, particularly through the lens of local covariance matrices.

6.1. Data Generating Process

We simulated $n=500$ observations with three covariates, $X_1,X_2,$ and $X_3$, independently drawn from a uniform distribution on $[0,1]$. The response variable was generated according to the nonlinear regression function

$$Y_i=m\left({\mathbf{X}}_i\right)+{\epsilon }_i,m\left({\mathbf{X}}_i\right)=sin\left(\pi X_{i1}{X}_{i2}\right)+cos\left(\pi X_{i3}\right),$$

where ${\epsilon }_i\sim N(0,0.1^2)$ represents independent additive noise. Although the population covariance matrix of $\mathbf{X}$ is approximately diagonal, the local covariance matrices ${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)$ vary across the covariate space due to finite sampling and kernel weighting. This variability is central to local polynomial fitting and directly impacts estimator stability and precision.

6.2. Missing Data Mechanisms

Missingness was introduced in both the response Y and the covariate $X_1$ using three distinct mechanisms. Under missing completely at random (MCAR), values of Y and $X_1$ were independently omitted with probabilities of 0.2, ensuring that missingness is independent of both observed and unobserved variables. Under missing at random (MAR), the probability of missingness in Y depended on the observed covariate $X_1$ via a logistic function, $P(R_Y=0\mid X_1)={\mathrm{logit}}^{-1}(0.5X_1-0.25)$, while $X_1$ remained MCAR with a 10% missing rate. Finally, under missing not at random (NMAR), the probability of missingness in Y depended on its own (unobserved) value through $P(R_Y=0\mid Y)={\mathrm{logit}}^{-1}(Y-0.5)$, introducing a systematic distortion in the observed covariance structure. These designs allow us to study the differential impact of missing data mechanisms on estimation and covariance geometry.

6.3. Estimation Procedure

Local linear regression ($q=1$) with a Gaussian kernel was employed to estimate $m\left(\mathbf{x}\right)$ across the covariate space. Two approaches were compared: complete-case (CC) analysis, which uses only fully observed data, and inverse probability weighting (IPW), which weights observed responses by the inverse of the estimated probability of being observed. Bandwidth selection for each estimator was performed via leave-one-out cross-validation to balance bias and variance.

6.4. Performance Metrics

The performance of each estimator was assessed over $B=1000$ Monte Carlo replicates using four key metrics. Bias was quantified as the average deviation of the estimator from the true regression function across all observations. Root mean squared error (RMSE) measured overall estimation accuracy by incorporating both bias and variance. Coverage evaluated the proportion of times 95% confidence intervals contained the true $m\left(\mathbf{x}\right)$, reflecting the adequacy of variance estimation. Finally, the condition number of local covariance matrices ${\Sigma }_{\mathbf{X}}\left(\mathbf{x}\right)$ served as a diagnostic of numerical stability and identifiability in local fitting.

6.5. Numerical Results

Table 1 summarizes the results. Under MCAR, complete-case estimates were essentially unbiased, and RMSE remained moderate, reflecting the reduction in effective sample size but minimal distortion of covariance geometry. IPW performed similarly, as weighting has little effect when missingness is independent of the covariates. Under MAR, CC estimates exhibited a noticeable bias because systematic missingness in Y depended on $X_1$, while IPW successfully corrected most of this bias, although slight variance inflation remained. Under NMAR, both CC and IPW displayed bias, consistent with the theoretical expectation that missingness dependent on unobserved Y alters the observed covariance structure in ways that cannot be fully corrected by IPW alone. Notably, regions of the covariate space with greater missingness corresponded to higher condition numbers in the local covariance matrices, confirming the direct influence of covariance geometry on estimation stability.

6.6. Covariance Geometry Diagnostics

Figure 1 depicts a heatmap of the condition number of local covariance matrices across the covariate space under MAR and NMAR. Elevated condition numbers indicate near-singularity, corresponding to regions where missingness has reduced directional variation and local polynomial fitting is prone to instability. These visualizations empirically support the notion that missing data reshape the local covariance geometry, propagating their effects to estimator bias, variance, and coverage.

6.7. Interpretation and Discussion

The simulation results demonstrate that missing data mechanisms influence nonparametric estimation through the geometry of the covariate space. Under MCAR, missingness primarily reduces sample size without substantially altering local covariance structure. Under MAR, selective missingness introduces directional distortions in local covariance matrices, increasing estimator bias for complete-case analysis, but these distortions can be mitigated through weighting. Under NMAR, missingness depending on unobserved outcomes induces fundamental deviations in observed covariance geometry, which standard correction methods cannot fully address. Condition number analysis effectively identifies regions where estimation instability is likely, linking theoretical insights on covariance geometry to practical diagnostics. These findings reinforce the central thesis that variance–covariance structure is a latent geometric object that underpins the identifiability, stability, and efficiency of multiple nonparametric regression in the presence of missing data.

7. Discussion

The framework developed in this paper demonstrates that the variance–covariance matrix plays a central role in multiple nonparametric regression when data are incomplete. Unlike in classical parametric models, where the covariance primarily summarizes linear dependencies, in nonparametric regression the covariance matrix functions as a latent geometric object underlying local polynomial and kernel smoothing procedures. It governs identifiability, estimator stability, variance, and the validity of inference, particularly in the presence of missing data.

Local covariance matrices determine whether the regression surface can be uniquely identified in a neighborhood of a covariate point. If these matrices are singular or nearly singular—conditions that can arise due to collinearity or selective missingness—local polynomial and kernel estimators become unreliable. Furthermore, the conditioning of local covariance matrices directly affects numerical stability: poorly conditioned matrices amplify sampling noise and can result in oscillatory or unstable fits. Missing data exacerbate these issues by selectively reducing variation along certain directions in the covariate space.

The geometry of the local variance–covariance structure also influences estimator efficiency and adaptive bandwidth selection. Regions of low variability require broader smoothing to reduce variance, whereas regions with higher variation allow finer resolution. Missing data modify this geometry, inflating variance and complicating optimal bandwidth selection. Similarly, accurate statistical inference depends on covariance-adjusted asymptotic variances. Distortions caused by missingness necessitate careful estimation of standard errors and confidence intervals to maintain correct coverage.

Analysis of different missing data mechanisms illustrates their mechanism-specific impacts. Under missing completely at random (MCAR), the observed covariance matrix remains unbiased, though effective sample size is reduced, slightly inflating estimator variance. Under missing at random (MAR), conditioning on observed covariates introduces systematic distortions, which can bias complete-case estimators unless weighting or augmentation is applied. Missing not at random (NMAR) poses the greatest challenge, as the observed covariance may diverge fundamentally from the population structure, rendering standard nonparametric estimators biased and inconsistent.

Adopting a covariance-geometric perspective provides a unifying conceptual framework for understanding both the successes and failures of smoothing-based methods under missing data. It clarifies why regions of the covariate space with low local variance or ill-conditioned covariance matrices are particularly susceptible to instability and bias, and why conventional bandwidth selection or variance estimation procedures may fail without adjustments for missingness-induced distortion. Condition numbers of local covariance matrices serve as practical diagnostics, identifying regions where estimation is prone to numerical instability.

This perspective also suggests concrete directions for future research. Explicit asymptotic theory could characterize how missing data affect the distribution of local covariance matrices and associated estimators. Covariance-aware bandwidth selection procedures can be developed to adapt smoothing to local geometric properties while accounting for incomplete observations. Integration of doubly robust or augmented estimation approaches can mitigate bias under MAR and NMAR by explicitly modeling covariance distortion. Extensions to high-dimensional and functional data are also warranted, where local covariance matrices may be singular or infinite-dimensional, necessitating regularization or dimension reduction strategies.

By emphasizing the latent geometric role of variance–covariance structures, this framework bridges classical nonparametric smoothing theory with modern missing data methodology. It provides both theoretical insight and practical guidance for designing robust, efficient, and stable nonparametric regression procedures in the presence of incomplete data, ensuring reliable estimation and inference across diverse applied settings.

8. Conclusions

In this paper, we have highlighted the pivotal role of the variance–covariance structure in multiple nonparametric regression when data are incomplete. Unlike classical parametric regression, where covariance primarily summarizes linear dependencies, in nonparametric settings the covariance matrix functions as an implicit geometric backbone. It governs local identifiability, estimator stability, variance inflation, and the validity of statistical inference. The local structure of covariates directly affects the performance of kernel and local polynomial estimators, and missing data can distort this geometry in ways that depend on the mechanism of missingness.

Our theoretical analysis, supported by extensive simulation studies, illustrates how different missing data mechanisms impact the covariance structure and consequently the performance of nonparametric estimators. Under missing completely at random (MCAR), observed covariance matrices remain unbiased approximations of the population covariance; however, the reduction in effective sample size inflates estimator variance and slightly reduces precision. Under missing at random (MAR), the covariance structure is systematically altered conditional on observed covariates, introducing bias if unaccounted for. Methods such as inverse probability weighting or augmentation effectively mitigate these distortions, restoring near-unbiased estimation. In contrast, under missing not at random (NMAR), missingness depends on unobserved outcomes, potentially producing fundamental deviations between observed and population covariance matrices. Standard nonparametric procedures may fail under NMAR, highlighting the need for additional modeling assumptions or doubly robust estimation strategies.

Simulation results reinforce this conceptual framework. Regions of the covariate space where local covariance matrices are ill-conditioned correspond to higher estimator variance, greater bias, and reduced confidence interval coverage. This confirms the geometric interpretation of covariance: selective missingness diminishes variation along particular directions, which directly affects the stability and reliability of local smoothing procedures.

Recognizing covariance as a latent geometric object offers a principled foundation for methodological development. It enables the design of covariance-aware bandwidth selection strategies that adapt to local variability, the implementation of robust estimation procedures such as weighting, imputation, or doubly robust approaches that correct for covariance distortion, and the creation of diagnostic tools that monitor condition numbers or visualize local covariance to identify regions of potential instability.

Looking ahead, this covariance-geometric perspective can be extended to more complex settings, including high-dimensional covariates where local covariance matrices may be ill-conditioned or singular, functional or longitudinal data requiring infinite-dimensional generalizations of local covariance, and formal asymptotic theory that explicitly incorporates the geometry of covariance under MCAR, MAR, and NMAR mechanisms.

In summary, understanding and leveraging the variance–covariance structure provides a unifying framework connecting classical nonparametric smoothing theory with modern missing data methodology. It enables the development of estimation, inference, and diagnostic tools that are robust, efficient, and reliable, facilitating principled analysis of incomplete-data problems across diverse scientific domains.