Deep Spatial Quantile Regression Insights

• Deep spatial quantile regression is a framework that estimates conditional quantiles in spatially dependent data using advanced neural architectures.
• It combines classical spatial quantile methods with deep learning techniques, such as ConvLSTM layers and multi-task objectives, to address issues like quantile crossing.
• The approach is applied across diverse fields—from urban traffic forecasting to environmental analysis—to uncover tail behaviors and treatment effects beyond conditional means.

Deep spatial quantile regression is the study of conditional quantiles or full conditional quantile processes for spatially indexed or spatio-temporal responses when dependence across locations cannot be ignored and, in the narrow machine-learning sense, when nonlinear spatial features are learned by deep architectures. In the literature represented here, that narrow sense includes multi-output multi-quantile deep learning for spatio-temporal tensors and a semiparametric neural framework for spatial quantile treatment effects, while the broader technical background includes spatial autoregressive quantile models, copula-based joint quantile regression, smooth-density Bayesian models, and high-dimensional nonparametric spatial quantile processes (Rodrigues et al., 2018, Gong et al., 2 Sep 2025, Chen et al., 2019, Deb et al., 2024). The common objective is to move beyond conditional expectations and estimate how dispersion, tails, and uncertainty vary across space, time, and covariate profiles.

1. Quantile targets and spatial structure

In scalar-response form, quantile regression replaces the mean target $E(g\mid X)=X\theta$ with a conditional quantile function such as $Q_G(\tau\mid X_i)=X_i\theta_\tau$, estimated by minimizing an asymmetric absolute-loss criterion. This shift is substantive rather than cosmetic: it permits heterogeneous covariate effects across the lower, middle, and upper parts of the conditional distribution, and it is explicitly described as robust to heteroskedasticity, distribution-free, and useful for detecting heterogeneous tail effects (Cartone et al., 2019).

For multivariate spatial outputs, the target can be generalized through spatial quantiles. One formulation defines, for a random vector $Y\in\mathbb{R}^p$, the spatial quantile

$$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$

with $u\in B^{p-1}$, and then conditions this object on covariates to obtain $Q(u,x)\in\mathbb{R}^p$. In that setting, $u\approx 0$ corresponds to central quantiles or the spatial median, whereas $\|u\|$ near $1$ corresponds to extreme quantiles; the direction parameter is also used to avoid the classical quantile crossing issue in a multivariate setting (Deb et al., 2024).

Once observations are spatially or spatio-temporally indexed, the conditional quantile is no longer purely a function of covariates in the usual iid sense. Neighboring regions may interact through knowledge spillovers, policy spillovers, labor market interactions, omitted regional effects, or spatial clustering, and spatio-temporal arrays may exhibit both recurrent temporal dependence and local spatial coupling. Deep spatial quantile regression therefore sits at the intersection of distributional regression and spatial modeling: its defining problem is not merely estimating $Q_Y(\tau\mid X)$, but estimating it when the conditioning information or the latent structure is itself spatially organized.

2. Spatial dependence in quantile regression before deep architectures

A central classical formulation is the spatial autoregressive quantile model

$Q_G(\tau\mid X_i)=X_i\theta_\tau$0

which augments quantile regression by conditioning on a spatial lag $Q_G(\tau\mid X_i)=X_i\theta_\tau$1. In the European regional growth application, this framework is used on 187 European NUTS2 regions over 1981–2009 to study conditional $Q_G(\tau\mid X_i)=X_i\theta_\tau$2-convergence in GDP per worker. The methodological argument is that standard quantile regression ignores spatial dependence and can be misspecified and potentially inconsistent when spatial lag dependence is present; because $Q_G(\tau\mid X_i)=X_i\theta_\tau$3 is endogenous, the paper adopts IVQR with instruments $Q_G(\tau\mid X_i)=X_i\theta_\tau$4. Empirically, conditional $Q_G(\tau\mid X_i)=X_i\theta_\tau$5-convergence is supported, convergence is faster in the lower tail than in the center or upper tail, and the spatial autoregressive coefficient $Q_G(\tau\mid X_i)=X_i\theta_\tau$6 is smaller in the lower part of the distribution and increases substantially in the upper tail, indicating stronger spatial dependence among better-performing regions (Cartone et al., 2019).

A related econometric line studies spatial quantile autoregressive regression simultaneously over several quantiles. In that setting,

$Q_G(\tau\mid X_i)=X_i\theta_\tau$7

and endogeneity of the spatial lag is handled by instrumental variables with $Q_G(\tau\mid X_i)=X_i\theta_\tau$8. Interquantile commonality is then encoded through fused adaptive LASSO or fused adaptive sup-norm penalties on adjacent quantile differences. The paper proves oracle properties and reports that the penalized procedures have smaller MedSE than conventional IV quantile regression, with FAL or FL often performing best and revealing piecewise-constant coefficient paths across quantiles (Dong et al., 2021).

Another route is fully generative joint quantile regression for spatial data. Here the model writes

$Q_G(\tau\mid X_i)=X_i\theta_\tau$9

with spatial dependence introduced through a copula process on the latent uniforms $Y\in\mathbb{R}^p$0. Gaussian and $Y\in\mathbb{R}^p$1 copula processes are used, the latter to capture tail dependence, and the model supports spatial quantile smoothing at new locations through the conditional copula of $Y\in\mathbb{R}^p$2. The framework distinguishes itself from fully spatially varying coefficient models: it retains a global quantile regression relationship while adjusting inference and prediction for spatial correlation. In PM$Y\in\mathbb{R}^p$3 data from 339 stations in the northeastern United States, the spatial joint models attained the smallest held-out check loss across quantiles, and WAIC strongly favored the spatial specifications over nonspatial JQR (Chen et al., 2019).

These pre-deep formulations established the main technical themes later inherited by neural models: quantile noncrossing, explicit spatial dependence, tail behavior, endogeneity, and the distinction between marginal quantile modeling and full generative modeling.

3. Deep spatio-temporal joint mean and quantile regression

A direct deep-learning contribution is DeepJMQR, a multi-output multi-quantile deep learning model for spatio-temporal tensors $Y\in\mathbb{R}^p$4. Its objective is to predict, for each grid cell $Y\in\mathbb{R}^p$5 and future time $Y\in\mathbb{R}^p$6, both the conditional mean $Y\in\mathbb{R}^p$7 and multiple conditional quantiles $Y\in\mathbb{R}^p$8. The architecture uses a stack of ConvLSTM layers so that temporal dependencies are handled by recurrent dynamics and spatial dependencies by convolutions in the recurrence; the hidden and cell states are tensors $Y\in\mathbb{R}^p$9. From the shared latent state $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$0, a single head predicts $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$1 outputs through hard parameter sharing,

$$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$2

equivalently as a $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$3 convolution over the latent tensor (Rodrigues et al., 2018).

Training uses a multi-task objective that combines squared error for the mean with the tilted loss or pinball loss for each quantile. The paper explicitly interprets this either as a sum of task losses or as a regularized mean-prediction objective in which quantile losses act as extra supervision; optimization is performed with automatic differentiation and Adam. Quantile crossing is not prevented by explicit monotonicity constraints. Instead, the model relies on multi-task coupling and hard parameter sharing, and the paper reports that this strategy yields zero or near-zero crossings while simultaneously significantly outperforming state-of-the-art quantile regression methods (Rodrigues et al., 2018).

The empirical evaluation spans a heteroscedastic motorcycle benchmark, NYC taxi demand, and Copenhagen/Nørrecampus traffic speeds. For NYC taxi demand, the data comprise 1.1 billion taxi trips in Manhattan discretized into a $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$4 grid with 30-minute bins and a 1-hour-ahead forecasting task. For Copenhagen traffic, the data comprise 6 months of crowdsourced traffic speed data on 9 consecutive road segments with 5-minute bins. Across experiments, DeepJMQR dramatically reduces crossings relative to independently trained quantile models; on Copenhagen, independent deep models show about 81k crossings while DeepJMQR has zero, and on NYC taxi the independent deep baseline has huge crossing loss and about 292k crossings while DeepJMQR reduces this to near zero. The paper also reports statistically significant improvements in MAE and RMSE for mean prediction on NYC and Copenhagen, intervals that are often narrower at similar coverage, and only a negligible increase in training time, with Copenhagen training taking about 28.19 minutes for DeepJMQR versus about 27.69 minutes for the mean-only model (Rodrigues et al., 2018).

4. Causal deep spatial quantile regression

A more explicitly spatial neural formulation appears in causal spatial quantile regression. The observational data are $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$5, where $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$6 is the outcome, $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$7 the treatment, $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$8 observed covariates, and $$Q(u)=\arg\min_{q\in\mathbb{R}^p} E\!\left[\|Y-q\|+\langle u,Y-q\rangle\right],$$9 spatial coordinates. The conditional density and CDF are modeled through spline mixtures,

$u\in B^{p-1}$0

with second-order $u\in B^{p-1}$1-spline and $u\in B^{p-1}$2-spline bases and softmax neural-network weights $u\in B^{p-1}$3. Spatial structure is enriched by multi-resolution basis functions $u\in B^{p-1}$4 built from a compactly supported Wendland basis, used in a manner described as similar to DeepKriging-style spatial deep learning (Gong et al., 2 Sep 2025).

The central estimand is the spatial quantile treatment effect

$u\in B^{p-1}$5

with a spatially averaged version obtained by averaging over $u\in B^{p-1}$6. Identification assumes SUTVA, consistency, ignorability, and positivity. The model is semiparametric: the output layer is constrained by spline basis functions, while feature extraction from $u\in B^{p-1}$7 is performed by a flexible neural network. This allows nonlinear covariate effects, treatment-covariate interactions, spatially varying effects, distributional asymmetry, and quantile-specific behavior within a single framework (Gong et al., 2 Sep 2025).

The paper also proposes an adjustment for spatial hidden confounders. At a target site $u\in B^{p-1}$8, one selects observations within a fixed neighborhood radius, fits the model locally on that subregion, and uses the estimated coefficients to predict counterfactuals across the full region. The stated rationale is that hidden confounders that are smoother than the treatment in space are approximately constant within a local neighborhood, so local fitting can reduce bias. Simulations are performed on $u\in B^{p-1}$9 spatial locations with $Q(u,x)\in\mathbb{R}^p$0 observations per location and $Q(u,x)\in\mathbb{R}^p$1 Monte Carlo replicates, with RMISE as the main metric. In the hidden-confounding scenario, Model 5 performs best, nearly as well as in the unconfounded case; spatial confounding adjustment helps Models 1 and 2 substantially, whereas when spatial features are already rich, local adjustment adds little (Gong et al., 2 Sep 2025).

The real-data application uses North Carolina birth records from 1988–2002 for first-time white mothers at 691 distinct ZIP-code centroid locations. Maternal smoking is estimated to have consistently negative effects across all birth-weight quantiles, with particularly severe impacts in the lower quantile regions; at $Q(u,x)\in\mathbb{R}^p$2, the smoking effect is negative throughout North Carolina but spatially heterogeneous in magnitude, with stronger or more visible heterogeneity in southern and northeastern regions (Gong et al., 2 Sep 2025).

5. High-dimensional, nonparametric, and fully generative extensions

Several adjacent methods are not deep-learning architectures but are central to the technical ecology of deep spatial quantile regression. One nonparametric spatio-temporal framework considers a process $Q(u,x)\in\mathbb{R}^p$3 observed over discrete times and a growing set of spatial locations, with response vector $Q(u,x)\in\mathbb{R}^p$4 and conditional spatial quantile process $Q(u,x)\in\mathbb{R}^p$5. The estimator minimizes a kernel-weighted multivariate spatial quantile criterion and is computed by an iteratively reweighted least squares scheme with weights $Q(u,x)\in\mathbb{R}^p$6; under $Q(u,x)\in\mathbb{R}^p$7, it attains

$Q(u,x)\in\mathbb{R}^p$8

and the paper develops asymptotic normality, simultaneous confidence bands, and Gumbel-limit tests based on a functional central limit theorem for martingale differences. In smart-meter electricity data from the Thames Valley Vision project, quantile curves are reported as smooth, with no visible quantile crossing, while tests of homogeneity for clusters 2 and 12 yield p-values from $Q(u,x)\in\mathbb{R}^p$9 to $u\approx 0$0 (Deb et al., 2024).

For image- and function-valued responses, high-dimensional spatial quantile function-on-scalar regression models

$u\approx 0$1

with coefficient functions in an RKHS and a Student-$u\approx 0$2 copula with Matérn spatial correlation. Estimation proceeds through penalized quantile regression and a primal-dual algorithm designed for large images, while the copula layer recovers the full conditional spatial distribution and supports conditional image generation. The paper proves minimax lower and upper bounds under fixed and random designs and reports applications to ADNI corpus callosum DTI data and hippocampal surface data, where age and disease effects vary by both spatial location and quantile level (Zhang et al., 2020).

A Bayesian smooth-density spatial quantile regression model takes a different route by parameterizing the entire conditional quantile function with I-spline basis functions in the center of the distribution and generalized Pareto tails in the extremes. The regression effects $u\approx 0$3 vary jointly with quantile level and space, monotonicity is enforced through positivity constraints on spline coefficients, and the induced density is guaranteed to be continuous or $u\approx 0$4-times differentiable under explicit threshold-matching constraints. Because the coefficient surfaces have Gaussian process priors, the model also induces predictor-dependent non-stationary covariance functions. In simulations it performs especially well under heavy tails, and in benzene measurements around a Corpus Christi refinery it shows that the 95th percentile effect of one source is especially strong and positive on northern fence-line sites (Brantley et al., 2019).

Other semiparametric precursors emphasize local smoothing and robustness rather than full generative structure. Functional-coefficient spatial quantile regression models

$u\approx 0$5

and estimates the unknown $u\approx 0$6 by local linear M-estimation under weak spatial mixing, covering both stationary and nonstationary fields with spatial trends. In the soil250 application, effects of soil chemistry variables and neighboring observations differ across the 0.15, 0.5, and 0.85 quantiles, going beyond mean or median-only analyses (Lu et al., 2014). For multivariate responses, local constant and local bilinear multiple-output quantile/depth regression estimate conditional halfspace depth contours by kernel weighting around a target covariate value; the local bilinear method is explicitly designed to adapt to nonlinear and heteroskedastic dependence and to recover the conditional contour system that characterizes the full conditional law (Hallin et al., 2015).

6. Empirical patterns and methodological distinctions

The application domains of this literature are unusually broad: European regional economic convergence, PM$u\approx 0$7 concentrations, refinery benzene measurements, soil chemistry, taxi demand, traffic speed, smart-meter electricity demand, neuroimaging, and birth weight all appear as case studies. The recurring empirical benefit is that spatial quantile models reveal effects that are weak or invisible at the conditional mean. In European growth, convergence is supported at all quantiles but is strongest in the lower tail and accompanied by stronger high-quantile spatial autocorrelation; in PM$u\approx 0$8, URB is more negative at lower and middle quantiles while L25E10 becomes pronounced at high quantiles; in the benzene study, upper-tail source effects are spatially heterogeneous; and in North Carolina birth data, maternal smoking is harmful across all quantiles but especially in the lower tail (Cartone et al., 2019, Chen et al., 2019, Brantley et al., 2019, Gong et al., 2 Sep 2025).

Several methodological distinctions structure the field. First, spatial dependence and spatially varying coefficients are not equivalent. Joint spatial quantile regression with latent copulas keeps a global quantile regression relationship and uses dependence in the latent $u\approx 0$9, whereas smooth-density and varying-coefficient models let coefficient functions themselves vary over space (Chen et al., 2019, Lu et al., 2014). Second, noncrossing is pursued by heterogeneous mechanisms: DeepJMQR relies on hard parameter sharing without explicit monotonicity constraints; joint Bayesian quantile models impose monotonicity analytically through quantile-function parameterization; multivariate spatial quantiles use a directional parameter $\|u\|$0 to avoid the classical quantile crossing issue in a multivariate setting (Rodrigues et al., 2018, Chen et al., 2019, Deb et al., 2024). Third, causal deep spatial quantile regression does not eliminate identification requirements: SUTVA, consistency, ignorability, and positivity remain explicit assumptions, and the treatment of spatial hidden confounders requires an additional local-neighborhood adjustment procedure (Gong et al., 2 Sep 2025).

A common misconception is to use “deep spatial quantile regression” as though it referred to a single methodological family. The papers considered here separate more sharply. Some are directly neural and representation-learning based; others are explicitly described as non-deep kernel, RKHS, Bayesian, or semiparametric frameworks that serve as precursors, contrasts, or statistical backbones for high-dimensional spatial quantile analysis (Deb et al., 2024, Zhang et al., 2020). This suggests that deep spatial quantile regression is best understood not as a replacement for classical spatial quantile statistics, but as a layer built on top of an already mature literature on spatial autoregression, copula processes, spline-based quantile functions, local M-estimation, and multivariate spatial quantiles.