Spatial processes have occupied center stage in statistical theory and applications for the last few decades. Their voracious use can largely be explained by geographically tagged data becoming increasingly commonplace in modern applications. Such data are often composed of complex variables which are no longer amenable to a Gaussian assumption. For example, spatially indexed counts [see for e.g., 61, 9, 2, 39], proportions [see, for e.g., 23, 26, 21, 20], time to event or, survival outcomes [see, for e.g. 4, 44, 67] are some frequently occurring variables where spatial processes have proved invaluable in performing uncertainty quantification. The purpose being to quantify unobserved dependence introduced within the variable of interest due to varying location. The cornerstone for such studies is a spatially indexed process variable of interest, often termed as a response process and denoted by, y ( s ). This is accompanied by covariate information X ( s ) = [ x 1 ( s ) , x 2 ( s ) , … , x p ( s ) ]. Here s ∈ S = { s 1 , s 2 , … , s L } is the spatial indexing and, S is a finite set of indices or, locations over which the response variable and covariates are observed. The investigator often encounters observation–level covariates that account for response specific characteristics when learning such processes. It becomes important to understand which of these covariates are important contributors to variation in the data. From a model parsimony standpoint, model choice becomes an important issue to the investigator. In statistical theory, this problem is often addressed by performing shrinkage or, variable selection on the model coefficients. Moreover, performing spatial uncertainty quantification produces accurate inference for model coefficients, which also raises the concern regarding a more “honest” subset of covariates within X ( s ) that primarily determine the variation in y ( s ). The crown jewel of our contribution is Bayesian methodology for performing spatial uncertainty quantification and model choice simultaneously.

Spatial process modeling generally requires a hierarchical specification of an unobserved random effect within the model [14]. Maintaining a hierarchical approach allows for exploration of the effect of covariates, X ( s ), and the random effect w ( s ) jointly on the response process y ( s ). Particularly, considering generalized linear spatial process modeling, it is assumed that y ( s ) ∣ β , w ( s ) arise in a conditionally independent fashion from a member of the exponential family with mean μ ( s ) such that g ( μ ( s ) ) = X ( s ) β + F ( s ) w ( s ), where g is a monotonic link function, β are model coefficients, or fixed effects and F ( s ) is a spatial incidence matrix. In contrast to Gaussian response processes where a direct hierarchical specification on the response is feasible, modeling a non-Gaussian spatial process leverages the generalized linear model framework to employ a latent process specification [see for e.g., 64, 17, 16, 65, 5]. This is facilitated by the existence of a valid joint probability distribution, π y ( s 1 ) , y ( s 2 ) , … , y ( s L ) ∣ β , θ p r , where θ p r denotes process parameters required for specification of w ( s ) [see discussion in section 6.2, 5]. This hierarchical specification gives us a natural way to perform variable selection by incorporating shrinkage priors into the hierarchical prior formulation. There are several choices of shrinkage priors that differ in their prior specification such as: the discrete spike-and-slab prior [see, for e.g, 46, 27], other priors based on the Gaussian-gamma family in linear Gaussian models [see, for e.g, 51, 7], the continuous spike-and-slab prior [see, for e.g, 33], the Bayesian counterparts of LASSO and elastic net [see, for e.g, 50, 41], mixtures of g–priors [see, for e.g, 42], the horseshoe priors and its variants [see, for e.g, 13], among several others. We use the continuous spike-and-slab prior to effect shrinkage on model coefficients.

We focus on a subset of probability distributions within the exponential family, termed as the exponential dispersion family [see, for e.g., 34, 35, 36, 37]. It allows the dispersion along with the mean to vary across observations, suppressing the need for having a constant dispersion across observations. We focus on a particular member of the family, the Tweedie compound Poisson-gamma (CP-g), more commonly referred to as the Tweedie (probability) distribution [see, 58]. The corresponding random variable is constructively defined as a Poisson sum of independently distributed Gamma random variables. Allowing for a varying dispersion across observations enables exploration of the effect of covariates X ( s ) on the mean and the dispersion separately, by employing two separate generalized linear models (GLMs). This gives rise to the double generalized linear model (DGLM) [see, for e.g., 54, 59, 56]. Hierarchical frameworks for specifications of DGLMs were first developed in [40]. Although not mandatory, it is customary to use the same covariates, X ( s ) for both the mean and dispersion GLMs to avoid ambiguities in model specification. Previous work [see for e.g., 29, 30] uses this approach and considers developing inference for DGLMs under a frequentist framework. Inference on spatial effects is obtained through penalizing the graph Laplacian. In this paper we adopt a Bayesian discourse by supplementing the DGLM framework with the continuous version of the spike-and-slab prior to effect shrinkage and thereby achieve better model specification. We integrate the spike and slab prior into our hierarchical prior formulation for both mean and dispersion models. We show that these priors provide a natural way of incorporating sparsity into the model, while offering straightforward posterior sampling in the context of our spatial DGLMs.

The scale for spatial indexing is assumed to be point-referenced. For example, latitude-longitude or easting–northing. Generally, s ∈ R 2 . Specification of a neighborhood structure or, proximity is naturally important when attempting to quantifying the behavior of response in locations that are near each other. We select the Euclidean distance between locations. This results in a Gaussian process [see, for e.g., 60] prior on the spatial process, w ( s ). Other choices exist for specifying such spatial process, w ( s ), for e.g. log-Gamma [see 10, 11] etc. We particularly focus on a Gaussian process specification on w ( s ) ∼ G P ( 0 , K ( · ) ), where K is a covariance function. For arbitrary locations, s and s ′ , dependence between y ( s ) and y ( s ′ ) is specified through K ( s , s ′ ), which governs the covariance between w ( s ) and w ( s ′ ). For point-referenced data, the Matérn family [see, for e.g., 45] provides the most generic and widely adopted covariance specification.

Next, we address Bayesian model specification. In the absence of such concerns for the hierarchical process models discussed above, prior specification follows the generic framework, [ data ∣ process , θ ˜ ] × [ process ∣ θ ˜ ] × [ θ m ] × [ θ p r ]. Here, θ ˜ = { θ m , θ p r } denote model parameters [see for e.g. 8, 24, 5, Chapter 6, p. 125]. In particular, θ p r constitute parameters instrumental in specification of the process, while θ m are other model parameters. We adopt a proper prior on θ p r to avoid the risk of generating improper posteriors [see, for e.g., 6]. Building a Bayesian variable selection framework that facilitates model specification for θ m requires an additional layer of hierarchical prior specification, appending the latter framework with variable selection parameters, θ v s and, thereby producing (1.1) [ data ∣ process , θ ] × [ process ∣ θ ] × [ θ p r ] × [ θ m ] × [ θ v s ] , where θ = { θ ˜ , θ v s }. We resort to Markov Chain Monte Carlo (MCMC) sampling [see, for e.g. 12, 28] for performing posterior inference on θ. The novelty of our approach lies in the simple Bayesian computation devised—employing only Gibbs sampling updates for θ v s . To the best of our knowledge hierarchical Bayesian frameworks for fitting (a) Tweedie DGLMs, (b) spatial Tweedie DGLMs with (or without) variable selection, do not exist in the statistical literature. Evidently, proposed methodology in (1.1) remedies that.

The ensuing developments in the paper are organized as follows: In Section 2 we detail the proposed statistical framework outlining Tweedie distributions—the likelihood and parameterization, model formulation and the hierarchical prior specification. Section 3 provides comprehensive synthetic experiments that document the efficacy of our proposed statistical framework for Bayesian variable selection in spatial DGLMs. Section 5 considers application of the developed framework to automobile insurance premiums in Connecticut, USA during 2008. Additional synthetic experiments capturing various performance aspects for the models are provided in the Supplementary Material.

The Tweedie distribution produces observations composed of exact zeros with a continuous Gamma tail. Their ability to model mixed data types featuring exact zeros and continuous measurements jointly makes them suitable for modeling response arising from a variety of domains. Some of the current applications include, actuarial science [see 55, 62, 29, 30], ecology [see for e.g., 57], public health [63], environment [see for e.g., 38], ecology [53], gene expression studies [43]. As evidenced by these applications, the presence of unobserved dependence between observations, affecting the quality of inference, is not unlikely. In the following subsections we provide more details on Tweedie distributions, followed by the model formulation and hierarchical prior specification.

We begin with an observation—the spatial heterogeneity that our models aim to quantify is not observed in real life. Hence, it is imperative to document the accuracy of estimating such effects through synthetic experiments. Settings used are outlined—we consider varying proportion of zeros (15%, 30%, 60%, 80% and 95%) under which the quality of posterior inference for θ is assessed. Proportion of zeros can be interpreted as an inverse signal-to-noise ratio for the synthetic response. For the sake of brevity, we only show the results for synthetic experiments pertaining to Bayesian variable selection in the presence of spatial effects. Additional simulations can be found in the online Supplement. To construct the synthetic data we consider three scenarios pertaining to model structure, (a) there is no overlap (i.e. selected β’s and γ’s do not intersect) (b) 50% overlap (in the union of all selected variables across the mean and dispersion models) (c) 100% overlap between mean and dispersion model specification. We use 10 covariates including an intercept, where the columns of the synthetic design matrices X and Z are hierarchically centered and scaled, independently sampled Gaussian variables with mean 0 and variance 1. Naturally, specification of the true β, w and γ parameters determine the proportion of zeros in the synthetic response. Table 2 contains the parameter specifications used. The true value of the index parameter, ξ = 1.5. In an attempt to produce a synthetic setup that resembles reality we simulate, w ∼ N ( 5 ( sin ( 3 π s 1 ) + cos ( 3 π s 2 ) ) , 1 ) (see Figure 1, second row). The alternative route would be to fix values of σ 2 and ϕ s and generate a realization w ∼ N L ( 0 L , σ 2 R ( ϕ s ) ) (see Figure 1, first row). Under each setting we consider M = 10 replications. Within each replication we fit all the proposed modeling frameworks as shown in Table 1. The hyper-parameter settings used while specifying priors for the models are, a ξ = 1, b ξ = 2, a σ = a σ β = a σ γ = 2, b σ = b σ β = b σ γ = 1 (producing inverse-Gamma priors with mean 1 and infinite variance), a ϕ s = 0, b ϕ s = 30, σ β − 2 = σ γ − 2 = 10 − 6 and ν 0 = 5 × 10 − 4 , producing a vague and non-informative hierarchical prior. We maintain an FDR of 5% for all settings while performing model selection. The sample size varies from N = { 2 × 10 3 , 5 × 10 3 , 1 × 10 4 } and the number of locations are L = 1 × 10 2 . Across replications, false positive rate (FPR) and true positive rate (TPR) are computed to measure accuracy of our model selection procedure. To record the quality of estimation, we compute the mean squared error (MSE), for e.g. M S E ( β ) = 1 p ∑ i β = 1 p ( β i β − β ˆ i β ) 2 , which can be computed similarly for the other parameters. We also compute average coverage probabilities, for e.g. considering these probabilities for β we define C P ( β ) = 1 M ∑ m = 1 M I ( β t r u e ∈ ( l m ( β ) , u m ( β ) ) ), where l m ( β ) and u m ( β ) are the lower and upper 95% HPD intervals respectively for β in replication m; we obtain coverage probabilities for w and γ similarly. The results obtained under the above settings are shown in Table 3. The first column is named configuration (abbreviated as config.) with entries denoting the proportion of overlap between selected coefficients in the mean and dispersion models, which is indicative of model structure. This is estimated by observing the overlap between selected variables following the model fit (for models M2 and M4). No variable selection is performed for models M1 and M3.

From the results shown, we see that models M1 and M2 perform poorly. Estimates, β ˆ remain fairly unaffected as compared to γ ˆ and ξ ˆ, where all of the variation not quantified, yet present in the synthetic data spills over to corrupt and compromise the quality of estimates. This also does not produce reliable results pertaining to model structure recovery for M1 and M2. However, significant improvements show up with M3 and M4. Particularly, under higher proportion of zeros in the synthetic data (low signal to noise ratio) the performance of M4 remains stable with respect to model structure recovery and estimation of parameters (refer to Table 3), thereby producing robust inference among the models in comparison. As an example within our simulation setting, under 95% of 0s in the data and under low sample sizes, for example 2000 or, 5000, the estimates of model coefficients and spatial effects in M3 and M4 are adversely affected by locations having fewer non-zero observations. This observation addresses the concern around specifying DGLMs without spatial random effects in a scenario where the data displays spatial variation. The results demonstrate expected gains when our model in its full capacity is used instead of an usual DGLM.

We use the MCMC algorithm featuring MALA updates for β , w and γ. Chain lengths are set to 1 × 10 4 , with the initial 5,000 samples as burn-in and thin the rest by selecting every 10-th sample which reduces any remaining auto-correlation and produces 500 independent posterior samples for each setting. The posterior estimate, θ ˆ is then obtained using the produced samples by computing the median or a MAP estimate as applicable for the model. Coverage probabilities for model M4 remained sufficiently high ( ≈ 1) across all settings; only declining marginally for w (remaining above 90 %) under high proportions of zeros (low signal to noise ratio) in the data. We performed additional synthetic experiments to showcase (a) the performance of M3 with respect to the quality of estimation for spatial effects and (b) the performance of M2. They are detailed in the Supplementary materials—we briefly outline its contents in the next section.

The online Supplement to this paper contains details on the derivations of the posteriors essential for constructing MCMC subroutines. They are outlined in Section S1. Section S2 features additional simulation experiments that supplement those outlined previously in Section 3. It documents performance of M2, shown in Table S1, contains results of experiments for scenarios featuring spatial covariates, shown in Table S2, and varying spatial patterns as seen in Figure 1, shown in Tables S3 and S4. Convergence diagnostics are shown for selected model parameters (index parameter ξ) in Section S2.2. Contents of the R-package are described in Section S2.3. Finally, results for models M1 and M3 pertaining to the real data analysis described in the next section appear in Section S3, Tables S5, S6, S7 and S8.

We analyze automobile insurance premiums for collision coverage in the state of Connecticut, for the year 2008. The data is obtained from a comprehensive repository named Highway Loss Data Institute (HLDI) maintained by the independent non-profit, Insurance Institute for Highway Safety (IIHS) (https://www.iihs.org/) working towards reducing losses arising from motor vehicle collisions in North America. We briefly describe the variables contained in HLDI data. It records covariate information at three levels for an insured individual. They are as follows, a.

Individual Level: (i) accident and model year of the vehicle, (ii) age, gender, marital status.

We denote payments towards a policy made by individual i, residing in zip-code j as y i j with an exposure of t i j . We assume that the policy-premium defined as y i j ∗ = y i j t i j ∼ T w μ i j , ϕ i j , ξ , which implies y i j ∼ T w ( t i j μ i j , t i j 2 − ξ ϕ i j , ξ ) using the scale invariance property. The following hierarchical DGLM is then specified, (5.1) log μ i j ( s i ) = − log t i j + x i j T ( s i ) β + f i j ( s i ) T w ( s i ) , log ϕ i j = − ( 2 − ξ ) log t i j + z i j T γ , when considered with a spatial specification, where the terms − log t i j and − ( 2 − ξ ) log t i j act as offsets for the respective mean and dispersion models. Given the covariates described at the beginning p = q = 29, producing a N t r × ( p − 1 ) design matrix for the mean model and a N t r × q design matrix for the dispersion model. The model in (5.1) specifies model M3 from Table 1, M4 is obtained by specifying π ( θ v s ) from (2.5) on θ m . M1 is obtained by setting f i j ( s i ) = 0 and M2 is obtained by specifying π ( θ v s ) on θ m for the resulting model. For M1 and M2 we include an intercept in the mean model. We fit models M1–4 on the training data. Model selection is performed using FDR based variable selection on the posterior MCMC samples obtained from fitting models M2 and M4, controlling for FDR at 1%. The performance of M1–4 is assessed using the Akaike Information Criteria (AIC) [3]. While specifying (2.4) and (2.5) the hyper-parameter settings used are, σ β 2 = σ γ 2 = 10 6 , a σ β = a σ = a σ γ = 2 and b σ β = b σ = b σ γ = 1 (generating inverse-gamma priors with mean 1 and infinite variance), a ϕ s = 0, b ϕ s = 60. We maintain a ξ = 1 and b ξ = 2 for all models. For ease of implementation, the fractal parameter, ν is fixed at 0.5, producing the exponential covariance kernel. We consider 1 × 10 5 MCMC iterations for generating samples from respective posteriors with burn-in diagnosed at 5 × 10 4 and include every 20-th sample to compute posterior estimates as our thinning strategy. Convergence is assessed through inspecting trace plots. Proposal variances were scaled in an adaptive fashion to provide optimal acceptance rates of 58% (MALA) and 33% (MH). Predictive performance for models M1–4 is judged based on square root deviance on the out-sample data. The results are shown in Table 6. Optimal values are marked in bold. We show the results for estimated model coefficients featuring Bayesian variable selection (models M2 and M4) in Tables 4 and 5. Results for M1 and M3 are postponed to Tables S5, S6, S7 and S8 in the Supplementary Materials. Posterior estimates for the spatial effects in models M3 and M4 are showed in Figure 4. Zip-codes with significant effects are color coded appropriately.

Comparing the models we observe that M4 produces the most optimal model fit criteria among the models considered. This extends to out-sample performance when predicting policy premiums. Plots produced for spatial effects in models M3 and M4 are mean adjusted. Since specification of M3 and M4 differ only in presence/absence of the hierarchical Bayesian variable selection component, the produced spatial effects mimic each other after adjusting for the mean. Comparing results in Tables 4 and 5 we observe that including the spatial effect results in more categories for vehicle age, driver age and deductible being selected. Overall, we observe that the findings remain consistent with our earlier research [see, 29] but within a more robust proposed model choice framework. This is evident when comparing model estimates between Table 5 and Table S7 and S8 in the Supplement. The marital status and interaction between gender and marital status is not selected. We conclude by observing that the spatial effects are significantly positive in major cities in Connecticut, indicating a higher spatial risk, as opposed to sparsely populated regions showing significantly lower risk.

Double generalized linear models have not seen much use after their inception by [40]. Hindrances presented by ambiguities existing around model specification/choice have been addressed in this paper. We propose Bayesian modeling frameworks that perform model selection using continuous spike and slab priors for hierarchical double generalized linear Tweedie spatial process models. Leveraging Langevin dynamics we are able to successfully produce practical implementations for the proposed frameworks which would otherwise remain unachievable with standard MCMC techniques. The proposed algorithms are available as a publicly accessible package for the R–statistical environment. Although the formulation considers the CP-g densities, evidently such modeling could be effected under any probabilistic framework that allows for varying dispersion. The application offers some key insights into the actuarial domain. It is generally believed that marital status and gender play a key role. However, the model inference suggests otherwise, with marital status not being selected as a significant feature.

Future work is aimed at extending this framework in multiple directions. Firstly, with the advent of modern Bayesian variable selection priors—for example, the Bayesian Lasso, the Horseshoe prior etc., a comparative model selection performance remains to be seen when considered within hierarchical DGLM formulations. Secondly, with the emerging techniques for handling large spatial and spatiotemporal data [see, for e.g., 31] the DGLM framework could be extended to model spatially or, spatio-temporally indexed observations over massive geographical domains. With respect to our application, this would allow us to investigate properties of the premium surface over much larger domains, for instance a country-wide study. Finally, extending these models to a spatiotemporal setting could be achieved using spatiotemporal covariance kernels that are commonly used. Depending on the nature of spatial and temporal interaction, we can have separable and non-separable kernels at our disposal [see, 15, and references therein]. Bayesian variable selection could then be effected to examine resulting changes in model specification upon inclusion of random effects that address spatiotemporal variation in the data.