# GPBoost: How Spatial Observations & Boosting Work Together Based on the Medium article by Fabio Sigrist (June 21, 2023) and the GPBoost source code analysis, here's the full explanation: --- ## 1. High-Level: How the Model Solves Spatial Observations GPBoost combines **two complementary components** into a single mixed-effects model: ### The Model Structure \[ y_i = \underbrace{f(\mathbf{x}_i)}_{\text{Tree-boosting (fixed effects)}} + \underbrace{u(s_i) + v_{cl(i)}}_{\text{Random effects (spatial + grouped)}} + \varepsilon_i \] - **\(y_i\)** — log(GDP/capita) for region \(i\) - **\(f(\mathbf{x}_i)\)** — Non-linear function of predictors (L, K, Pop, Edu), learned via tree-boosting - **\(u(s_i)\)** — Spatial random effect at location \(s_i\), modeled by a Gaussian Process with an exponential covariance function - **\(v_{cl(i)}\)** — Grouped random effect for the cluster (there are 2 clusters) - **\(\varepsilon_i\)** — i.i.d. Gaussian noise ### How it handles spatial structure 1. **Each region** \(i\) has centroid coordinates (Longitude, Latitude) 2. A Gaussian Process prior is placed on the spatial random effects: \(\mathbf{u} \sim \mathcal{N}(0, \Sigma)\) 3. The covariance matrix \(\Sigma\) captures spatial correlation — regions that are geographically close have more similar random effects 4. The grouped random effect captures additional structure (two clusters of regions that may share unobserved common factors) 5. Both the GP covariance parameters (\(\sigma^2, \rho, \tau^2\)) and the grouped random effect variance are estimated via **maximum likelihood** or **REML** This means the model automatically accounts for spatial autocorrelation beyond what the predictors can explain. --- ## 2. What Observations Are Mapped Into the Covariance Function The covariance function operates on **the 242 region centroids** — each defined by (Longitude, Latitude) coordinates. Specifically: - **`gp_coords = data[, c("Long", "Lat")]`** — the centroid spatial coordinates of every European region - The exponential covariance between two regions \(i\) and \(j\) is: \[ \text{Cov}(u(s_i), u(s_j)) = \sigma^2 \cdot \exp\left(-\frac{d(s_i, s_j)}{\rho}\right) + \tau^2 \cdot \mathbf{1}_{i=j} \] Where: - \(d(s_i, s_j)\) = Euclidean distance between centroid coordinates of regions \(i\) and \(j\) - \(\sigma^2\) = marginal variance (sill) — estimated by ML - \(\rho\) = range parameter — controls how quickly correlation decays with distance, estimated by ML - \(\tau^2\) = nugget — measurement error / micro-scale variation on the diagonal, estimated by ML **Crucially**, the 242 spatial locations produce a **242 × 242 covariance matrix** (not 484 × 484) — because each region appears twice in the dataset (years 2000 and 2021), but the spatial random effect \(u(s_i)\) is the same for a given region across time. The **grouped random effect** \(v_{cl(i)}\) and the **boosting fixed effects** handle the within-region temporal variation. Additionally, **grouped random effects** are mapped via `group_data = data[, c("cl")]` — the cluster ID variable that assigns each region to one of 2 clusters. --- ## 3. How Boosting Comes Into Play ### The Key Innovation of GPBoost The article's opening line says it directly: > *"The GPBoost algorithm extends linear mixed effects and Gaussian process models by **replacing the linear fixed effects function** with a **non-parametric non-linear function modeled using tree-boosting**."* In a traditional linear mixed model: \[ y = \underbrace{X\beta}_{\text{linear fixed effects}} + Zu + \varepsilon \] GPBoost replaces \(X\beta\) with \(f(X)\), a **gradient boosted tree ensemble**: \[ y = \underbrace{f(X)}_{\text{tree-boosting}} + Zu + \varepsilon \] ### How it works in the demo - **Predictors** used in the tree-ensemble: L (employment share), K (fixed capital), Pop (population), Edu (tertiary education) - **Boosting parameters**: `learning_rate = 0.01`, `max_depth = 2`, `num_leaves = 2^10`, `min_data_in_leaf = 10`, `lambda_l2 = 0`, `nrounds = 37` - The tree ensemble is trained using gradient boosting — each tree fits the **partial residuals** (what's left after accounting for the random effects) ### The Training Algorithm The training alternates between two estimation tasks: 1. **Boosting step**: Given current estimates of the random effects \(\mathbf{u}\) and \(\mathbf{v}\), the tree ensemble \(f(X)\) is grown by gradient boosting to minimize the log-likelihood (or equivalently, to fit the response adjusted for random effects) 2. **Random effects (GP) step**: Given the current tree ensemble \(f(X)\), the Gaussian process parameters (\(\sigma^2, \rho, \tau^2\)) and the grouped random effect variance are re-estimated via maximum likelihood, and the random effects are predicted (via best linear unbiased prediction / BLUP) This is essentially a **coordinated optimization** where boosting handles the non-linear signal from the predictors, and the GP handles the spatial correlation in the residuals. ### Why combine them? - **Pure tree-boosting** ignores spatial dependence → spatial autocorrelation in residuals violates independence assumptions and hurts generalization - **Pure spatial GP model** with linear fixed effects → misses non-linear relationships between predictors and the response - **GPBoost** → gets the best of both: non-linear prediction via boosting + proper spatial uncertainty quantification via Gaussian processes --- ## Summary Table | Component | What it does | What observations map into it | |---|---|---| | **Exponential covariance function** | Models spatial correlation between regions based on geographic distance | 242 region centroids (Long, Lat) — computes Euclidean distances between all pairs | | **Grouped random effects** | Captures cluster-level shared variation | Cluster ID `cl` (2 clusters) — each region assigned to a cluster | | **Tree-boosting ensemble** | Learns non-linear relationships between predictors and GDP | 484 rows × 4 predictors (L, K, Pop, Edu) — all data points | | **GP parameter estimation** | Covariance parameters (\(\sigma^2, \rho, \tau^2\)) estimated via ML | Uses all 484 data points through the likelihood function |