Spatiotemporal Autoregressive Gamma Process

• The spatiotemporal autoregressive gamma process is defined as a Bayesian model for count data using latent gamma frailties with explicit stationarity constraints.
• It introduces latent Poisson decompositions and sparse neighborhood dependencies to achieve computational efficiency with linear scaling in large datasets.
• Empirical evaluations demonstrate accurate parameter estimation and superior predictive performance compared to traditional CAR and INLA-based models.

The spatiotemporal autoregressive gamma process is a Bayesian modeling strategy for spatiotemporal count data in which observed counts are conditionally Poisson and the latent positive spatiotemporal frailties evolve through a gamma-process autoregression that is guaranteed stationary across the time dimension under explicit parameter restrictions. In the formulation of Cheng and Li, the construction replaces the more common latent Gaussian specification for log-rates with a conjugate Gamma–Poisson hierarchy, introduces latent Poisson decompositions, and yields efficient posterior sampling with sparse spatial dependence and effectively linear computational scaling in the total number of space–time observations (Cheng et al., 26 Jul 2025).

1. Formal definition and hierarchical construction

Let $s_1,\dots,s_m$ denote $m$ spatial locations, $t=1,\dots,T$ time periods, and $x_t(s_i)$ a $p\times 1$ covariate at $(s_i,t)$. The observed counts are modeled as

$$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$

where $U_t(s_i)>0$ is the spatiotemporal frailty.

The defining feature is the conjugate Gamma–Poisson hierarchy placed on $U_t(\cdot)$ so as to induce an autoregressive gamma process in time. For $t=2,\dots,T$ and $m$0, the latent decomposition is

$m$1

and, for each of the $m$2 neighbors $m$3,

$m$4

Conditionally,

$m$5

while the initial layer is

$m$6

An equivalent one-line representation is

$m$7

with

$m$8

This is described as the non-central AR(1) Gamma process in time.

This construction is motivated by the limitations of previous models that decompose logarithms of the response Poisson rates into fixed effects and spatial random effects, where the latter is typically assumed to follow a latent Gaussian process, the conditional autoregressive model, or the intrinsic conditional autoregressive model. Since log-Gaussian is not conjugate to Poisson, those implementations must resort to approximation methods like INLA or Metropolis moves on latent states in MCMC algorithms and exhibit several approximation and posterior sampling challenges. The autoregressive gamma specification is designed precisely to avoid that non-conjugacy.

2. Stationarity, spatial dependence, and sparsity

The model imposes a local spatial dependence structure. For each $m$9, one chooses $t=1,\dots,T$0 of size $t=1,\dots,T$1, for example the $t=1,\dots,T$2 nearest neighbors, and prespecifies nonnegative weights $t=1,\dots,T$3 on $t=1,\dots,T$4 summing to $t=1,\dots,T$5. Within this structure, the $t=1,\dots,T$6 term models own-location temporal autocorrelation, and the $t=1,\dots,T$7 term smooths toward neighbors.

The temporal stationarity condition is explicit. If $t=1,\dots,T$8 and $t=1,\dots,T$9, then for any nonnegative neighbor-weight matrix $x_t(s_i)$0 whose rows sum to $x_t(s_i)$1, the AR–Gamma process in $x_t(s_i)$2 is stationary (Cheng et al., 26 Jul 2025). This stationarity guarantee is central because it is built into the process definition rather than imposed only indirectly through a transformed latent Gaussian layer.

Sparsity is equally central. Because each row of the weight matrix $x_t(s_i)$3 has at most $x_t(s_i)$4 nonzero entries, consisting of self plus $x_t(s_i)$5 neighbors, all conditionals become sparse and require $x_t(s_i)$6 per full sweep. Since $x_t(s_i)$7, this sparse design is what makes large-$x_t(s_i)$8 spatiotemporal count analysis computationally feasible within a fully Bayesian posterior-sampling framework.

A common misconception in this area is that strong spatiotemporal dependence for counts must be mediated through Gaussian random effects on the log-scale. The present formulation provides a direct counterexample: dependence is carried by positive frailties under a Gamma–Poisson hierarchy, with temporal persistence, neighborhood smoothing, and stationarity all expressed in the latent gamma process itself.

3. Latent decomposition and posterior computation

The key computational device is the introduction of latent Poisson variables $x_t(s_i)$9, with $p\times 1$0 for the self-link and $p\times 1$1 for neighbor links. Under this augmentation, the full conditionals become standard.

For $p\times 1$2, the frailty update is

$p\times 1$3

Similar Gamma updates apply at $p\times 1$4 and $p\times 1$5, omitting the next or previous $p\times 1$6 term.

For each $p\times 1$7, $p\times 1$8, and $p\times 1$9,

$(s_i,t)$0

where $(s_i,t)$1 and $(s_i,t)$2 depend on $(s_i,t)$3. In practice one uses the Devroye rejection sampler for the Poisson–Bessel mixture.

The remaining global-parameter updates are

$(s_i,t)$4

$(s_i,t)$5

and

$(s_i,t)$6

The regression coefficient update is the sole non-conjugate component. Its posterior is

$(s_i,t)$7

with $(s_i,t)$8. No closed form is available, so the algorithm uses a Metropolis–Hastings random-walk with Hessian-based adaptive proposal.

The resulting Gibbs sampler proceeds by sampling all $(s_i,t)$9 in $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$0, all $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$1 in $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$2, then $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$3 in $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$4, and finally the Metropolis step for $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$5 in $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$6. Because $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$7 and $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$8, the total cost is effectively $$Y_t(s_i)\mid U_t(s_i),\beta \sim \operatorname{Pois}\{U_t(s_i)\cdot \exp[x_t(s_i)'\beta]\},$$9.

4. Predictive distribution and extrapolation in space and time

Posterior prediction is handled by composition sampling. For prediction at future times $U_t(s_i)>0$0 and/or new spatial sites $U_t(s_i)>0$1, the latent frailties are propagated via

$U_t(s_i)>0$2

with

$U_t(s_i)>0$3

The counts are then generated from

$U_t(s_i)>0$4

This predictive mechanism is operationally important because it applies both to temporal forecasting and to prediction at new spatial locations. The model description states that it delivers satisfactory performance in predicting at new spatial locations and time intervals (Cheng et al., 26 Jul 2025). A plausible implication is that the same latent augmentation that supports efficient posterior sampling also makes posterior predictive simulation straightforward, since the transition kernel remains in the gamma/noncentral-gamma family.

5. Empirical behavior and comparative performance

The reported simulation study is extensive. For $U_t(s_i)>0$5 up to $U_t(s_i)>0$6, $U_t(s_i)>0$7 up to $U_t(s_i)>0$8, and over a wide range of $U_t(s_i)>0$9, the posterior means of $U_t(\cdot)$0 track the true values with negligible bias, and the Mean Absolute Error (MAE) of fitted $U_t(\cdot)$1 remains low (Cheng et al., 26 Jul 2025).

Relative to six MCMC models in CARBayesST and two INLA-based models, the AR–Gamma model is reported to have comparable or lower MAE in-sample and out-of-sample, similar or better information criteria (DIC, WAIC), vastly higher effective sample size per CPU second, full conjugacy except $U_t(\cdot)$2 (one Metropolis), and linear scaling in $U_t(\cdot)$3.

Real-data analyses are reported for weekly COVID-19 cases and deaths with approximately $U_t(\cdot)$4 countries and $U_t(\cdot)$5 weeks. The recorded performance is MAE approximately $U_t(\cdot)$6–$U_t(\cdot)$7 for cases and approximately $U_t(\cdot)$8–$U_t(\cdot)$9 for deaths, MAPE approximately $t=2,\dots,T$0–$t=2,\dots,T$1, together with clear superiority over CAR and INLA in cross-validation.

These findings are significant primarily because the comparison class includes both MCMC-based CAR models and INLA-based approximations. The reported results therefore address two distinct questions at once: posterior efficiency and predictive adequacy. The model is presented as obtaining satisfactory model fitting, accurate parameter estimation, and strong out-of-sample prediction without abandoning full Bayesian posterior computation.

6. Relation to autoregressive gamma-process research

The spatiotemporal autoregressive gamma process for lattice-indexed count data is closely related to a broader gamma-process literature. The paper explicitly references the autoregressive gamma process for time series of counts due to Creal, Chib, and Shephard. In the spatiotemporal construction, that time-series logic is extended by combining self-location persistence, neighbor-weighted dependence, and sparse spatial coupling.

A related arXiv development is the measure-valued autoregressive gamma process of Bassetti, Casarin, and Iacopini, used as the latent intensity in a spatiotemporal shot-noise Cox process (Bassetti et al., 2023). There, the latent process is a Markov chain of gamma random measures on a Polish space $t=2,\dots,T$2, with one-step transition

$t=2,\dots,T$3

and an AR(1)+innovation representation

$t=2,\dots,T$4

If $t=2,\dots,T$5 are constant and $t=2,\dots,T$6, then $t=2,\dots,T$7, and the limiting law is invariant.

The two constructions differ in both inferential target and computational machinery. The spatiotemporal count model of Cheng and Li works directly with counts $t=2,\dots,T$8 and latent frailties $t=2,\dots,T$9, exploits latent Poisson decomposition, and uses a Gibbs sampler with a Metropolis step only for $m$00. By contrast, the shot-noise Cox process uses latent gamma random measures $m$01, a kernelized intensity

$m$02

and Bayesian inference via Particle-Gibbs with blocking and conditional Sequential Monte Carlo (Bassetti et al., 2023).

This suggests a useful conceptual distinction. The lattice-based spatiotemporal autoregressive gamma process is tailored to large spatiotemporal count arrays with sparse neighborhood structure and linear-in-$m$03 computation, whereas the measure-valued formulation is designed for point-process intensities, product densities, and pair-correlation analysis. Both, however, rely on noncentral-gamma transitions and on explicit stationarity conditions rather than latent Gaussian approximations.

7. Interpretation, scope, and methodological significance

Within Bayesian spatiotemporal modeling, the main significance of the spatiotemporal autoregressive gamma process lies in the conjunction of three features: conjugacy with Poisson observations, stationarity across time under transparent constraints, and sparse spatial dependence with effectively linear complexity. The model is therefore not merely a replacement likelihood or a computational trick; it is a distinct latent-process specification for positive spatiotemporal frailties.

Its scope is specifically spatiotemporal count data. The observed process is Poisson conditional on latent frailties and regression effects, and the latent dependence is encoded through the parameters $m$04, $m$05, the neighbor sets $m$06, and the weights $m$07. Because prediction at future times and new spatial sites is built from posterior samples of the same transition law, the inferential and predictive components are tightly aligned.

An important methodological point is that the model is not fully conjugate in every parameter: $m$08 still requires a Metropolis–Hastings random-walk with Hessian-based adaptive proposal. Accordingly, claims of efficiency rest not on universal closed-form updating, but on the fact that all full conditionals are standard except for $m$09, together with sparse conditionals and the latent Poisson decomposition. That distinction matters when comparing the model to latent Gaussian alternatives or to measure-valued gamma-process Cox-process models.

Taken together, the available results position the spatiotemporal autoregressive gamma process as a gamma-driven alternative to log-Gaussian, CAR, and ICAR-based spatiotemporal count models, with direct relevance for large count arrays, posterior prediction at new sites and times, and scalable Bayesian computation (Cheng et al., 26 Jul 2025).