Spatial Proportional Hazards Model with Differential Regularization

The Proportional Hazards (PH) model is one of the most widely used models in survival analysis, typically assuming a log-linear relationship between covariates and the hazard function. However, in the context of spatial survival data, where the time-to-event variable is associated with a spatial location within a given domain, this assumption is often unrealistic in capturing spatial effects. Thus, this paper proposes modeling the location effect through a nonparametric function of spatial location. The function is approximated using finite element methods on a triangulated mesh to accommodate irregular domains. Estimation is carried out within the classical partial likelihood framework, with smoothness of the spatial effect enforced through differential penalization. Using sieve methods, we establish the consistency and asymptotic normality of the parametric component. Simulations and two empirical applications demonstrate superior performance compared to existing approaches.

💡 Research Summary

The paper addresses a fundamental limitation of the classical Cox proportional hazards model when applied to spatial survival data: the assumption of a log‑linear covariate effect fails to capture complex spatial variation in the hazard. To overcome this, the authors introduce a non‑parametric spatial term h(p) that depends on the two‑dimensional location p of each observation. The hazard is written as
λ(t | x, p) = λ₀(t) exp{xᵀβ + h(p)}.
The spatial function h is assumed to belong to the Sobolev space H²(Ω) with homogeneous Neumann boundary conditions and a centering constraint ∫Ωh(p)dp = 0, which guarantees identifiability up to an additive constant.

Estimation proceeds within the partial‑likelihood framework, preserving the Cox model’s key advantage of treating the baseline hazard λ₀(t) as a nuisance parameter. To enforce smoothness of h, a differential penalty λ∫Ω(Δh)²dp is added, where Δ denotes the Laplacian. This “second‑order” penalty penalizes curvature rather than just roughness, leading to a more physically plausible spatial surface.

Because h is infinite‑dimensional, the authors adopt a sieve approach: they construct a sequence of finite‑dimensional subspaces that become dense in H²(Ω) as the sample size grows. The subspaces are built using finite element methods (FEM) on a triangulated mesh of the study region Ω. Importantly, they employ C¹‑conforming elements (e.g., Argyris elements) that guarantee the required second‑order differentiability of the basis functions, allowing the Laplacian‑squared penalty to be evaluated exactly. Each mesh is characterized by a size parameter η; the corresponding FEM space consists of piecewise quintic polynomials with 21 degrees of freedom per triangle, ensuring H²‑conformity.

The penalized partial‑likelihood objective is
Qₙ(β, h) = (1/n)∑_{i=1}ⁿ δ_i