1. Introduction

Reliable estimates of above-ground biomass density (AGBD) are essential for forest carbon accounting, REDD+, and conservation policies. Although Mexico's National Forest Inventory (INFyS) provides the most reliable field-based AGBD data, its limited spatial coverage and high acquisition costs constrain large-scale monitoring, particularly in tropical forests. Earth Observation (EO) products, such as ESA CCI Biomass and GEDI-derived canopy height, offer continuous spatial coverage but are affected by non-stationary relationships with AGBD and sensor saturation. Consequently, integrating EO and field data requires models that explicitly account for spatial dependence and uncertainty, making spatially varying coefficient (SVC) geostatistical models a suitable approach.

2. Research Gap
Models fitted over large, heterogeneous domains often borrow strength across broad ecological gradients and tend to shrink local predictions toward the national mean. While statistically efficient at the national scale, this can lead to regional miscalibration, particularly where Earth Observation (EO) covariates lose sensitivity. Both ESA Climate Change Initiative Biomass (CCI Biomass), which becomes weakly informative beyond roughly 150–200 Mg ha$^{-1}$, and GEDI-derived canopy height, which saturates near the upper end of tall tropical canopies, degrade precisely in dense tropical forests. Quintana Roo, with semi-evergreen tropical forest, secondary vegetation, coastal gradients, karstic substrates, and heterogeneous land-use histories, exemplifies this condition. A state-scale application therefore requires an independent modeling framework that jointly addresses response scale, data quality, and local covariate behavior, explicitly identifying where each EO source drives predictions.
3. Methodology

We develop such a reformulation for Quintana Roo within the SVC–SPDE framework. On a transformed response, the model is written in the standard additive form

$$ g{Y(\mathbf{s})}=[\alpha+\tilde{\alpha}(\mathbf{s})]+[\beta+\tilde{\beta}(\mathbf{s})]x_1(\mathbf{s})+[\eta+\tilde{\eta}(\mathbf{s})]x_2(\mathbf{s})+\epsilon(\mathbf{s}), $$

where $\tilde{\alpha},\tilde{\beta},\tilde{\eta}$ are mean-zero Matérn Gaussian fields approximated by the SPDE over a triangular mesh. The database combines INFyS plot AGBD, CCI Biomass, and GLAD-GEDI canopy height, harmonized around 2019, with covariates extracted as footprint-weighted means per plot. Because spatial covariance is inferred from the data, a single high-leverage record can distort local coefficient surfaces, so we applied a quality-control protocol—coordinate and unit verification, removal of duplicates and implausible values, and restriction to mainland forest. The SPDE mesh was confined to Quintana Roo rather than a national grid, improving computational efficiency.

The choice of $g$ matters once predictions return to the Mg ha$^{-1}$ scale, motivating our departure from the square-root response used nationally. Writing $\mu(\mathbf{s})$ for the linear predictor, the square-root link inverts by squaring, whereas the logarithmic link inverts by exponentiation:

$$ g=\sqrt{\cdot}: Y=\mu^2=\dots+2,\beta(\mathbf{s})\eta(\mathbf{s}),x_1x_2+\dots, \qquad g=\log: Y=e^{\alpha(\mathbf{s})},e^{\beta(\mathbf{s})x_1},e^{\eta(\mathbf{s})x_2}. $$

Squaring introduces covariate interaction terms, so the marginal effect of each covariate depends on the other and on the intercept, making EO effects entangled. The logarithm instead factorizes biomass into independent multiplicative terms, matching its power-law structure and preserving per-covariate effects. Zero-valued records were dropped, as $\log 0$ is undefined; since $\exp{E[\log Y]}\neq E[Y]$, predictions were back-transformed from posterior predictive samples. Accuracy was assessed with $R^2$, RMSE, MAE, and bias, calibration with leave-one-out Probability Integral Transform, and prediction structure through covariate-contribution surfaces $C_{\mathrm{CCI}}(\mathbf{s})=\beta(\mathbf{s})x_1(\mathbf{s})$ and $C_{\mathrm{GEDI}}(\mathbf{s})=\eta(\mathbf{s})x_2(\mathbf{s})$, with relative dominance $D_{\mathrm{CCI}}=|C_{\mathrm{CCI}}|/(|C_{\mathrm{CCI}}|+|C_{\mathrm{GEDI}}|)$.

4. Findings

The logarithmic formulation substantially outperformed the square-root one, returning an $R^2$ of 0.68 against 0.38, RMSE of 26.7 against 46.9 Mg ha$^{-1}$, MAE of 21.9 against 33.2 Mg ha$^{-1}$, and bias of 4.2 against 16.8 Mg ha$^{-1}$—a 78% gain in explained variance, a 43% reduction in RMSE, and a near four-fold reduction in relative bias, from 27.6% to 7.0%. These gains reflect the positive, right-skewed, heteroscedastic nature of biomass: aligning the linear predictor with multiplicative mechanisms prevents spatially varying coefficients from absorbing distortions arising from an unsuitable response scale. Applied locally, the PIT diagnostic confirms better-calibrated behaviour once variance is stabilised and anomalous records are removed.

Because the log scale keeps covariate effects separable, the contribution surfaces admit a clean reading that the squared scale does not. A coefficient map alone is misleading—a steep local slope contributes little where its covariate is near zero—so weighting each slope by the value its covariate takes reveals which EO source governs prediction at a given point. Mapping dominance over Quintana Roo exposes a spatial division of labor: zones where CCI Biomass carries the signal, where canopy height does, and where neither does and the latent field absorbs the residual baseline. This matters because similar canopy heights can correspond to different biomass densities depending on composition and degradation. Hunka et al. (2024) map spatially varying coefficients but stop there; carrying the analysis through to contribution and dominance turns a surface of parameters into an operational map of per-sensor reliability—and is coherent only under a response scale that does not entangle the covariates

5. Contributions

This study proposes an independent log-scale Spatially Varying Coefficient–Stochastic Partial Differential Equation (SVC–SPDE) formulation for tropical biomass that integrates quality control, a state-restricted mesh, and posterior-sample back-transformation into a unified workflow, improving performance in high-biomass regimes where Earth Observation (EO) covariates saturate. Its main contribution is the spatial decomposition of covariate effects using a logarithmic response that preserves separability on the natural scale, enabling identification of where each EO source governs prediction. The framework also supports posterior predictive aggregation over administrative and conservation units with explicit uncertainty, providing a transferable approach for biomass estimation.