# Bayesian Variable Selection in Double Generalized Linear Tweedie Spatial
Process Models

Aritra Haldera,∗†, Shariq Mohammedb,c† and Dipak K. Deyd

aDepartment of Biostatistics, Dornsife School of Public Health,
Drexel University, Philadelphia, PA, USA
bDepartment of Biostatistics, Boston University School of Public Health,
Boston University, Boston, MA, USA
cRafik B. Hariri Institute for Computing and Computational Science &
Engineering,
Boston University, Boston, MA, USA
dDepartment of Statistics, University of Connecticut, Storrs, CT, USA
∗corresponding author. E-mail<EMAIL_ADDRESS>
†equal contribution

###### Abstract

Double generalized linear models provide a flexible framework for modeling
data by allowing the mean and the dispersion to vary across observations.
Common members of the exponential dispersion family including Gaussian,
compound Poisson-gamma, Gamma, inverse-Gaussian are known to admit such
models. However, the lack of their use can be attributed to ambiguities that
exist in model specification under a large number of covariates and,
complications that arise when data from a chosen application display complex
spatial dependence. In this work we consider a hierarchical specification for
these models with a spatial random effect. The spatial effect is targeted at
performing uncertainty quantification by modeling dependence within the data
arising from location based indexing of the response. We focus on a Gaussian
process specification for the spatial effect. Simultaneously, we tackle the
problem of model specification for such hierarchical spatial process models
using Bayesian variable selection. It is effected through a continuous spike
and slab prior placed on the model parameters (or fixed effects). The novelty
of our contribution lies in the Bayesian frameworks developed for such models,
which have not been explored previously. We perform various synthetic
experiments to showcase the accuracy of our frameworks. These developed
frameworks are then applied to analyse automobile insurance premiums in
Connecticut.

Keywords: Bayesian Modeling, Gaussian Process, Hierarchical Spatial Models,
Spike and Slab Priors, Tweedie Double Generalized Linear Models.

##  1 Introduction

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., Wolpert and Ickstadt,, 1998; Best et al.,, 2000;
Agarwal et al.,, 2002; Lawson,, 2018), proportions (see for e.g., Gelfand,,
2000; Gelfand et al.,, 2005; Finley et al.,, 2009; Eidsvik et al.,, 2012),
time to event or, survival outcomes (see for e.g. Banerjee and Carlin,, 2004;
Martino et al.,, 2011; Zhou and Hanson,, 2015) 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(\mbox{\boldmath{$\mathbf{s}$}})$. This is accompanied by covariate
information
$\mbox{\boldmath{$\mathbf{X}$}}(\mbox{\boldmath{$\mathbf{s}$}})=[\mbox{\boldmath{$\mathbf{x}$}}_{1}(\mbox{\boldmath{$\mathbf{s}$}}),\mbox{\boldmath{$\mathbf{x}$}}_{2}(\mbox{\boldmath{$\mathbf{s}$}}),\ldots,\mbox{\boldmath{$\mathbf{x}$}}_{p}(\mbox{\boldmath{$\mathbf{s}$}})]$.
Here $\mbox{\boldmath{$\mathbf{s}$}}\in{\cal
S}=\\{\mbox{\boldmath{$\mathbf{s}$}}_{1},\mbox{\boldmath{$\mathbf{s}$}}_{2},\ldots,\mbox{\boldmath{$\mathbf{s}$}}_{L}\\}$
is the spatial indexing and, ${\cal 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
$\mbox{\boldmath{$\mathbf{X}$}}(\mbox{\boldmath{$\mathbf{s}$}})$ that
primarily determine the variation in $y(\mbox{\boldmath{$\mathbf{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 (Clark and Gelfand,, 2006).
Maintaining a hierarchical approach allows for exploration of the effect of
covariates, $\mbox{\boldmath{$\mathbf{X}$}}(\mbox{\boldmath{$\mathbf{s}$}})$
and, the random effect
$\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}})$ jointly on
the response process $y(\mbox{\boldmath{$\mathbf{s}$}})$. Particularly,
considering generalized linear spatial process modeling, it is assumed that
$y(\mbox{\boldmath{$\mathbf{s}$}})\mid\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}})$
arise in a conditionally independent fashion from a member of the exponential
family with mean $\mu(\mbox{\boldmath{$\mathbf{s}$}})$ such that
$g(\mu(\mbox{\boldmath{$\mathbf{s}$}}))=\mbox{\boldmath{$\mathbf{X}$}}(\mbox{\boldmath{$\mathbf{s}$}})\mbox{\boldmath$\beta$}+\mbox{\boldmath{$\mathbf{F}$}}(\mbox{\boldmath{$\mathbf{s}$}})\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}})$,
where $g$ is a monotonic link function, $\beta$ are model coefficients, or
fixed effects and
$\mbox{\boldmath{$\mathbf{F}$}}(\mbox{\boldmath{$\mathbf{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., Zeger and Karim,, 1991;
Diggle et al.,, 1998; Dey et al.,, 2000; Zhang,, 2002; Banerjee et al.,,
2014). This is facilitated by the existence of a valid joint probability
distribution,
$\pi\left(y(\mbox{\boldmath{$\mathbf{s}$}}_{1}),y(\mbox{\boldmath{$\mathbf{s}$}}_{2}),\ldots,y(\mbox{\boldmath{$\mathbf{s}$}}_{L})\mid\mbox{\boldmath$\beta$},\mbox{\boldmath$\theta$}_{pr}\right)$,
where $\mbox{\boldmath$\theta$}_{pr}$ denotes process parameters required for
specification of
$\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}})$ (see
discussion in section 6.2, Banerjee et al.,, 2014). 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,
Mitchell and Beauchamp,, 1988; George and McCulloch,, 1993), other priors
based on the Gaussian-gamma family in linear Gaussian models (see, for e.g,
Raftery et al.,, 1997; Berger et al., 2001b, ), the continuous spike-and-slab
prior (see, for e.g, Ishwaran and Rao,, 2005), the Bayesian counterparts of
LASSO and elastic net (see, for e.g, Park and Casella,, 2008; Li and Lin,,
2010), mixtures of $g$–priors (see, for e.g, Liang et al.,, 2008), the
horseshoe priors and its variants (see, for e.g, Carvalho et al.,, 2010),
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.,
Jørgensen,, 1986, 1987, 1992; Jorgensen,, 1997). 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, Tweedie et al.,,
1984). 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 $\mbox{\boldmath{$\mathbf{X}$}}(\mbox{\boldmath{$\mathbf{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., Smyth,, 1989; Verbyla,, 1993; Smyth and Verbyla,,
1999). Hierarchical frameworks for specifications of DGLMs were first
developed in Lee and Nelder, (2006). Although not mandatory, it is customary
to use the same covariates,
$\mbox{\boldmath{$\mathbf{X}$}}(\mbox{\boldmath{$\mathbf{s}$}})$ for both the
mean and dispersion GLMs to avoid ambiguities in model specification. Previous
work (see for e.g., Halder et al.,, 2021, 2022) 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,
$\mbox{\boldmath{$\mathbf{s}$}}\in\mathbb{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., Williams and Rasmussen,, 2006) prior on the spatial
process, $\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}})$.
Other choices exist for specifying such spatial process,
$\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}})$, for e.g.
log-Gamma (see Bradley et al.,, 2018, 2020) etc. We particularly focus on a
Gaussian process specification on
$\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}})\sim
GP(\mbox{\boldmath{$\mathbf{0}$}},K(\cdot))$, where $K$ is a covariance
function. For arbitrary locations, $\mathbf{s}$ and
$\mbox{\boldmath{$\mathbf{s}$}}^{\prime}$, dependence between
$y(\mbox{\boldmath{$\mathbf{s}$}})$ and
$y(\mbox{\boldmath{$\mathbf{s}$}}^{\prime})$ is specified through
$K(\mbox{\boldmath{$\mathbf{s}$}},\mbox{\boldmath{$\mathbf{s}$}}^{\prime})$,
which governs the covariance between $w(\mbox{\boldmath{$\mathbf{s}$}})$ and
$w(\mbox{\boldmath{$\mathbf{s}$}}^{\prime})$. For point-referenced data, the
Matérn family (see, for e.g., Matérn,, 2013) 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,
$[\text{data}\mid\text{process},\widetilde{\mbox{\boldmath$\theta$}}]\times[\text{process}\mid\widetilde{\mbox{\boldmath$\theta$}}]\times[\mbox{\boldmath$\theta$}_{m}]\times[\mbox{\boldmath$\theta$}_{pr}]$.
Here, $\widetilde{\mbox{\boldmath$\theta$}}=\\{\mbox{\boldmath$\theta$}_{m}$,
$\mbox{\boldmath$\theta$}_{pr}\\}$ denote model parameters (see for e.g.
Berliner,, 2000; Gelfand and Banerjee,, 2017; Banerjee et al.,, 2014, Chapter
6, p. 125). In particular, $\mbox{\boldmath$\theta$}_{pr}$ constitute
parameters instrumental in specification of the process, while
$\mbox{\boldmath$\theta$}_{m}$ are other model parameters. We adopt a proper
prior on $\mbox{\boldmath$\theta$}_{pr}$ to avoid the risk of generating
improper posteriors (see, for e.g., Berger et al., 2001a, ). Building a
Bayesian variable selection framework that facilitates model specification for
$\mbox{\boldmath$\theta$}_{m}$ requires an additional layer of hierarchical
prior specification, appending the latter framework with variable selection
parameters, $\mbox{\boldmath$\theta$}_{vs}$ and, thereby producing

$[\text{data}\mid\text{process},\mbox{\boldmath$\theta$}]\times[\text{process}\mid\mbox{\boldmath$\theta$}]\times[\mbox{\boldmath$\theta$}_{pr}]\times[\mbox{\boldmath$\theta$}_{m}]\times[\mbox{\boldmath$\theta$}_{vs}],$
(1)

where
$\mbox{\boldmath$\theta$}=\\{\widetilde{\mbox{\boldmath$\theta$}},\mbox{\boldmath$\theta$}_{vs}\\}$.
We resort to Markov Chain Monte Carlo (MCMC) sampling (see, for e.g. Carlin
and Louis,, 2008; Girolami and Calderhead,, 2011) for performing posterior
inference on $\theta$. The novelty of our approach lies in the simple Bayesian
computation devised—employing only Gibbs sampling updates for
$\mbox{\boldmath$\theta$}_{vs}$. 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) 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.

##  2 Statistical Framework

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 Smyth and Jørgensen,, 2002; Yang et al.,,
2018; Halder et al.,, 2021, 2022), ecology (see for e.g., Swallow et al.,,
2016), public health (Ye et al.,, 2021), environment (see for e.g., Kokonendji
et al.,, 2021), ecology (Shono,, 2008), gene expression studies (Mallick et
al.,, 2022). 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.

###  2.1 The Exponential Dispersion Family: Tweedie Distributions

The Tweedie family of distributions belong to the exponential dispersion (ED)
family of models whose probability density/mass function has the generic form,

$\pi(y\mid\theta,\phi)=a(y,\phi)\exp\left\\{\phi^{-1}(y\theta-\kappa(\theta))\right\\},$
(2)

where $\theta$ is the natural or canonical parameter, $\phi$ is the dispersion
parameter, and $\kappa(\theta)$ is the cumulant function. Characterizing the
Tweedie family is an index parameter $\xi$, varying values of which produce
different members of the family. For e.g., the CP-g is obtained with
$\xi\in(1,2)$, for $\xi=1$ we obtain a Poisson and $\xi=2$ produces a Gamma
distribution, for $\xi\in(0,1)$ they do not exist (for further details see
Table 1, Halder et al.,, 2021). We are particularly interested in the CP-g
distributions in this work. In general, for the ED family we have the mean,
$\mu=E(y)=\kappa^{\prime}(\theta)$ and the variance,
$Var(y)=\phi\kappa^{\prime\prime}(\theta)$. For the CP-g we have
$\kappa(\theta)=(2-\xi)^{-1}\\{(1-\xi)\theta\\}^{2-\xi/1-\xi}$. Using the
relation, $\kappa^{\prime}(\theta)=\mu$, some straightforward algebra yields,
$\kappa(\theta)=(2-\xi)^{-1}\mu^{2-\xi}$ and
$\kappa^{\prime\prime}(\theta)=\mu^{\xi}$, implying $Var(y)=\phi\mu^{\xi}$ and
denoting $\alpha=(1-\xi)^{-1}(2-\xi)$ we have,

$a(y,\phi)=1\cdot
I(y=0)+y^{-1}\sum_{j=1}^{\infty}\left[\frac{y^{-\alpha}(\xi-1)^{\alpha}}{\phi^{1-\alpha}(2-\xi)}\right]^{j}\frac{1}{j!\Gamma(-j\alpha)}I(y>0).$

Evidently, $\pi(0\mid\theta,\phi)=\exp\\{-\phi^{-1}\kappa(\theta)\\}$.

We introduce some notation. Let $y_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})$
denote the $j$-th response at the $i$-th location
$\mbox{\boldmath{$\mathbf{s}$}}_{i}\in{\cal S}$, where $j=1,2,\ldots,n_{i}$
and $i=1,2,\ldots,L$ with $\sum_{i=1}^{L}n_{i}=N$. Together we denote,
$\mbox{\boldmath{$\mathbf{y}$}}=\mbox{\boldmath{$\mathbf{y}$}}(\mbox{\boldmath{$\mathbf{s}$}})=\\{\\{y_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\\}_{j=1}^{n_{i}}\\}_{i=1}^{L}$
as the $N\times 1$ response. Similarly,
$\mbox{\boldmath$\mu$}=\mbox{\boldmath$\mu$}(\mbox{\boldmath{$\mathbf{s}$}})=\\{\\{\mu_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\\}_{j=1}^{n_{i}}\\}_{i=1}^{L}$
and $\mbox{\boldmath$\phi$}=\\{\\{\phi_{ij}\\}_{j=1}^{n_{i}}\\}_{i=1}^{L}$
denotes the mean and dispersion vectors respectively. If
$\mbox{\boldmath{$\mathbf{y}$}}\mid\mbox{\boldmath$\mu$},\mbox{\boldmath$\phi$},\xi$
arises independently from a CP-g distribution, then the likelihood is given by

$\pi(\mbox{\boldmath{$\mathbf{y}$}}\mid\mbox{\boldmath$\mu$},\mbox{\boldmath$\phi$},\xi)=\prod_{i=1}^{L}\prod_{j=1}^{n_{i}}a_{ij}(y_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\mid\phi_{ij})\times\exp\left[\phi_{ij}^{-1}\left(\frac{y_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\mu_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})^{1-\xi}}{1-\xi}-\frac{\mu_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})^{2-\xi}}{2-\xi}\right)\right].$
(3)

Working with the likelihood, $\pi(\cdot)$, when devising computation,
evaluating the infinite series representation of $a(y,\phi)$ is required. The
two commonly used methods are—saddle-point approximation (see for e.g. Nelder
and Pregibon,, 1987; Smyth and Jørgensen,, 2002; Dunn and Smyth,, 2005;
Zhang,, 2013) and Fourier inversion (see, for e.g. Dunn and Smyth,, 2008). The
saddle-point approximation to (2) uses a deviance function based
representation where,
$\widetilde{\pi}(y\mid\mu,\phi)=b(y,\phi)\exp\left\\{-(2\phi)^{-1}d(y\mid\mu)\right\\}$.
For CP-g distributions, the deviance function is
$d(y\mid\mu)=d(y\mid\mu,\xi)=2\left\\{\left(y^{2-\xi}-y\mu^{1-\xi}\right)(1-\xi)^{-1}-\left(y^{2-\xi}-\mu^{2-\xi}\right)(2-\xi)^{-1}\right\\}$,
and $b(y\mid\phi,\xi)=(2\pi\phi y^{\xi})^{-1/2}I(y>0)+1\cdot I(y=0)\approx
a(y,\phi)\exp\left\\{\phi^{-1}y^{2-\xi}(1-\xi)^{-1}(2-\xi)^{-1}\right\\}$. We
performed experiments which showed that the saddle-point approximation
performs well when fewer zeros are present in the data. Under higher
proportions of zeros its performance was sub-optimal. However, albeit its
computationally intensive nature, in all scenarios the Fourier inversion based
method had stable performance. Hence, in this paper we use the evaluation of
$a(y,\phi)$ that is based on Fourier inversion. The adopted Bayesian approach
requires MCMC computation that relies on accurate likelihood evaluations.
Hence, we emphasize the importance of choosing the appropriate likelihood
function for application purposes. We denote the likelihood in (3) as
$Tw(\mbox{\boldmath$\mu$},\mbox{\boldmath$\phi$},\xi)$. Tweedie distributions
are the only members of the ED family that possess a scale invariance property
(see, for e.g., Jorgensen,, 1997, Theorem 4.1). This suggests for $c_{ij}>0$,
$\mbox{\boldmath{$\mathbf{y}$}}^{*}(\mbox{\boldmath{$\mathbf{s}$}})=\\{c_{ij}y_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\\}\sim
Tw(\mbox{\boldmath{$\mathbf{c}$}}^{T}\mbox{\boldmath$\mu$},{\mbox{\boldmath{$\mathbf{c}$}}^{2-\xi}}^{T}\mbox{\boldmath$\phi$},\xi)$
allowing observations with different scales of measurement to be modeled
jointly.

###  2.2 Model Formulation

Formulating DGLMs with spatial effects theoretically involves specification of
a spatial random effect in both mean and dispersion models. In such a scenario
complex dependencies can be specified to account for varied degrees of
uncertainty quantification. In the simplest case the corresponding spatial
random effects for the mean and dispersion models are independent Gaussian
processes. More complex scenarios can feature dependent Gaussian processes,
where the dependence arises from a cross-covariance matrix. Spatial random
effects in the mean model are readily interpretable—risk faced (adjustment to
mean premium paid in our case) owing to location. However, spatial random
effects in the dispersion model are not readily interpretable. Subsequently,
we choose to include spatial random effects only in the mean model for this
work.

Let $\mbox{\boldmath{$\mathbf{x}$}}_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})$
denote a $p\times 1$ vector of observed covariates for the mean model and
$\beta$ be the corresponding a $p\times 1$ vector of coefficients.
$\mbox{\boldmath{$\mathbf{f}$}}_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})$
denotes a $L\times 1$ vector specifying the location and
$w(\mbox{\boldmath{$\mathbf{s}$}}_{i})$ is the spatial effect at location
$\mbox{\boldmath{$\mathbf{s}$}}_{i}$ with
$\mbox{\boldmath{$\mathbf{w}$}}=\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}})=(w(\mbox{\boldmath{$\mathbf{s}$}}_{1}),w(\mbox{\boldmath{$\mathbf{s}$}}_{2}),\ldots,w(\mbox{\boldmath{$\mathbf{s}$}}_{L}))^{T}$
denoting the $L\times 1$ vector of spatial effects;
$\mbox{\boldmath{$\mathbf{z}$}}_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})$
denotes a $q\times 1$ vector of known covariates for the dispersion model and
$\gamma$ is the corresponding $q\times 1$ vector of coefficients. A Bayesian
hierarchical double generalized linear model (DGLM) using a non-canonical
logarithmic link function is specified as

$\log\mu_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})=\mbox{\boldmath{$\mathbf{x}$}}_{ij}^{T}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\mbox{\boldmath$\beta$}+\mbox{\boldmath{$\mathbf{f}$}}_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})^{T}\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}}_{i}),~{}\log\phi_{ij}=\mbox{\boldmath{$\mathbf{z}$}}_{ij}^{T}\mbox{\boldmath$\gamma$},$
(4)

which implies
$\mu_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})=\mu_{ij}(\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}})=\exp\left(\mbox{\boldmath{$\mathbf{x}$}}_{ij}^{T}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\mbox{\boldmath$\beta$}+\mbox{\boldmath{$\mathbf{f}$}}_{ij}^{T}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\right)$
and
$\phi_{ij}=\phi_{ij}(\mbox{\boldmath$\gamma$})=\exp(\mbox{\boldmath{$\mathbf{z}$}}_{ij}^{T}\mbox{\boldmath$\gamma$})$.

###  2.3 Hierarchical Prior Specification

In this section, we first present the hierarchical prior formulation for the
model and process parameters, $\tilde{\mbox{\boldmath$\theta$}}$, followed by
the prior formulation for the variable selection parameters
$\mbox{\boldmath$\theta$}_{vs}$. Prior specification for model and process
parameters are as follows:

$\begin{split}&\text{Model Parameters: }\xi\sim
U(a_{\xi},b_{\xi});\mbox{\boldmath$\beta$}\sim
N_{p}\left(\mbox{\boldmath{$\mathbf{0}$}}_{p},\mbox{\boldmath$\lambda$}_{\beta}^{T}\mbox{\boldmath{$\mathbf{I}$}}_{p}\right);\mbox{\boldmath$\gamma$}\sim
N_{q}\left(\mbox{\boldmath{$\mathbf{0}$}}_{q},\mbox{\boldmath$\lambda$}_{\gamma}^{T}\mbox{\boldmath{$\mathbf{I}$}}_{q}\right),\\\
&\text{Process Parameters: }\phi_{s}\sim
U\left(a_{\phi_{s}},b_{\phi_{s}}\right);\sigma^{-2}\sim{\rm
Gamma}\left(a_{\sigma},b_{\sigma}\right);~{}\nu\sim U(a_{\nu},b_{\nu}),\\\
&\text{Process: }\mbox{\boldmath{$\mathbf{w}$}}\sim
N_{L}\left(\mbox{\boldmath{$\mathbf{0}$}}_{L},\sigma^{2}\mbox{\boldmath{$\mathbf{R}$}}(\phi_{s})\right),\end{split}$
(5)

where $\mbox{\boldmath{$\mathbf{0}$}}_{m}$ is the $m\times 1$ zero vector and
$\mbox{\boldmath{$\mathbf{I}$}}_{m}$ is the $m\times m$ identity matrix;
$\mbox{\boldmath{$\mathbf{X}$}}\in\mathbb{R}^{n\times p}$ and
$\mbox{\boldmath{$\mathbf{Z}$}}\in\mathbb{R}^{n\times q}$ are design matrices
corresponding to the mean and dispersion models, with coefficients
$\mbox{\boldmath$\beta$}\in\mathbb{R}^{p}$ and
$\mbox{\boldmath$\gamma$}\in\mathbb{R}^{q}$, respectively;
$\mbox{\boldmath{$\mathbf{F}$}}\in\mathbb{R}^{n\times L}$ is the spatial
incidence matrix, $\mbox{\boldmath{$\mathbf{w}$}}\in\mathbb{R}^{L}$ is the
spatial effect, and
$\mbox{\boldmath{$\mathbf{R}$}}(\phi_{s})=\sigma^{2}(\phi||\Delta||)^{\nu}K_{\nu}(\phi||\Delta||)$,
where $K_{\nu}$ is the modified Bessel function of order $\nu$ (Abramowitz et
al.,, 1988), is the Matérn covariance kernel. Here
$\\{||\Delta||\\}_{ii^{\prime}}=||\mbox{\boldmath{$\mathbf{s}$}}_{i}-\mbox{\boldmath{$\mathbf{s}$}}_{i^{\prime}}||_{2}$,
the Euclidean distance between locations $\mbox{\boldmath{$\mathbf{s}$}}_{i}$
and $\mbox{\boldmath{$\mathbf{s}$}}_{i^{\prime}}$. $U(\cdot\mid a,b)$ denotes
the uniform distribution,
$N_{m}(\cdot\mid\mbox{\boldmath{$\mathbf{0}$}},\Sigma)$ is the $m$-dimensional
Gaussian with zero mean and covariance matrix $\Sigma$, and ${\rm
Gamma}(\cdot\mid a,b)$ is the Gamma distribution with shape-rate
parameterization. Note that the priors on
$\mbox{\boldmath$\lambda$}_{\beta}=(\lambda_{\beta,1},\ldots,\lambda_{\beta,p})$
and
$\mbox{\boldmath$\lambda$}_{\gamma}=(\lambda_{\gamma,1},\ldots,\lambda_{\gamma,q})$
are part of the variable selection priors and are discussed next. Referring to
the framework in (1), the resulting posterior from (5) establishes the
$[\text{data}\mid\text{process},\mbox{\boldmath$\theta$}]\times[\text{process}\mid\widetilde{\mbox{\boldmath$\theta$}}]$
step. Conditional posteriors for $\widetilde{\mbox{\boldmath$\theta$}}$ are
outlined in Section 7 of Supplementary Materials.

For the continuous spike-and-slab prior formulation,
$\mbox{\boldmath$\theta$}_{vs}=\\{\zeta_{\beta},\zeta_{\gamma},\sigma^{2}_{\beta},\sigma^{2}_{\gamma},\alpha_{\beta},\alpha_{\gamma}\\}$
(see, for e.g., Ishwaran and Rao,, 2005). Note that we have separate prior
formulations for mean and dispersion models. Let
$\mbox{\boldmath$\beta$}=(\beta_{1},\beta_{2},\ldots,\beta_{p})^{T}$ and
$\mbox{\boldmath$\gamma$}=(\gamma_{1},\gamma_{2},\ldots,\gamma_{q})^{T}$ be
the model coefficients corresponding to the mean and the dispersion models.
Let us define $\lambda_{\beta,u}=\zeta_{\beta,u}\sigma_{\beta,u}^{2}$ and
$\lambda_{\gamma,v}=\zeta_{\gamma,v}\sigma_{\gamma,v}^{2}$ for $u=1,\ldots,p$
and $v=1,\ldots,q$, respectively. We consider the following prior formulation:

$\displaystyle\pi(\mbox{\boldmath$\theta$}_{vs})=\begin{cases}\zeta_{\beta,u}\stackrel{{\scriptstyle
iid}}{{\sim}}(1-\alpha_{\beta})\delta_{\nu_{0}}(\cdot)+\alpha_{\beta}\delta_{1}(\cdot),~{}\zeta_{\gamma,v}\stackrel{{\scriptstyle
iid}}{{\sim}}(1-\alpha_{\gamma})\delta_{\nu_{0}}(\cdot)+\alpha_{\gamma}\delta_{1}(\cdot),\\\
\alpha_{\beta}\sim U(0,1),~{}\alpha_{\gamma}\sim U(0,1),\\\
\sigma_{\beta,u}^{-2}\stackrel{{\scriptstyle iid}}{{\sim}}{\rm
Gamma}(a_{\sigma_{\beta}},b_{\sigma_{\beta}}),~{}\sigma_{\gamma,v}^{-2}\stackrel{{\scriptstyle
iid}}{{\sim}}{\rm Gamma}(a_{\sigma_{\gamma}},b_{\sigma_{\gamma}}),\end{cases}$
(6)

where $\beta_{u}$ and $\gamma_{v}$ have normal priors with mean 0 and variance
$\zeta_{\beta,u}\sigma_{\beta,u}^{2}$ and
$\zeta_{\gamma,u}\sigma_{\gamma,u}^{2}$, respectively. Here,
$\delta_{c}(\cdot)$ denotes the discrete measure at $c$; hence,
$\zeta_{\beta,u}$ and $\zeta_{\gamma,u}$ are indicators taking values $1$ or
$v_{0}$ (small number close to 0) based on the selection of their
corresponding covariates. The probabilities of these indicators taking the
value 1 is given by $\alpha_{\beta}$ and $\alpha_{\gamma}$ respectively. We
place a uniform prior on these selection probabilities and an inverse-Gamma
prior on the parameters $\sigma_{\beta,u}^{2}$ and $\sigma_{\gamma,u}^{2}$.
The choice of the shape and rate parameters
($a_{\sigma_{\beta}},b_{\sigma_{\beta}};a_{\sigma_{\gamma}},b_{\sigma_{\gamma}}$)
of these inverse-Gamma priors induces a continuous bimodal distributions on
$\zeta_{\beta,u}\sigma_{\beta,u}^{2}$ and
$\zeta_{\gamma,u}\sigma_{\gamma,u}^{2}$ with a spike at $\nu_{0}$ and a right
continuous tail. Combining the priors in (5) and (6) completes the
hierarchical prior formulation for parameters $\theta$ as defined in (1).
Evidently, the above prior formulation allows for sufficient flexibility
regarding variations in implementation. For instance, a hierarchical Bayesian
framework for a simple DGLM can be obtained by omitting the process
specification and variable selection. Analogously, DGLMs featuring variable
selection or, spatial effects are obtained by omitting respective components
from the prior specification outlined previously.

Table 1: Proposed Bayesian Hierarchical Double Generalized Linear Modeling
Frameworks

Models | Frameworks | Specification ($\theta$) | Number of Parameters
---|---|---|---
M1 | DGLM | $\mbox{\boldmath$\theta$}_{m}$ | $p+q+1$
M2 | DGLM + Variable Selection | $\mbox{\boldmath$\theta$}_{m}$, $\mbox{\boldmath$\theta$}_{vs}$ | $3p+3q+1$
M3 | DGLM + Spatial Effect | $\mbox{\boldmath$\theta$}_{m}$, $\mbox{\boldmath$\theta$}_{pr}$ | $p+q+L+4$
M4 | DGLM + Spatial Effect + Variable Selection | $\mbox{\boldmath$\theta$}_{m}$, $\mbox{\boldmath$\theta$}_{vs}$, $\mbox{\boldmath$\theta$}_{pr}$ | $3p+3q+L+4$

###  2.4 Bayesian Estimation and Inference

In its full capacity (a model with spatial effects and variable selection) the
model structure with prior specifications in (5) and (6) contains $3p+L+3q+4$
parameters. Depending on the dimensions of
$\mbox{\boldmath{$\mathbf{X}$}}(\mbox{\boldmath{$\mathbf{s}$}})$ and
$\mathbf{Z}$, and the number of locations $L$, posterior inference can be a
sufficiently daunting task. Traditional Metropolis–Hasting (M-H) random walk
strategies are sub-optimal, involving costly pilot runs to determine viable
initial starting points and unreasonably long chains while performing MCMC
sampling. To avoid such issues, we use an adaptive rejection sampling while
leveraging the log-concavity of the posteriors to perform effective inference
that is not plagued by the above described issues (for more details, see, for
e.g, Girolami and Calderhead,, 2011). In the following, we describe (i)
briefly, our adaptive rejection MCMC sampling approach (more details are
provided in the Supplementary Materials), (ii) the identifiability issues on
the overall intercept that arise due to inclusion of a spatial effect and a
strategy to address this, and (iii) a false discovery rate (FDR)–based
approach for performing variable selection.

The joint posterior
$\pi(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath{$\mathbf{y}$}})$ generated as
a result of the hierarchical priors in eq. 5 is sampled using a hybrid
sampling strategy that includes M-H random walk and the Metropolis-Adjusted
Langevin Algorithm (MALA) (Roberts and Stramer,, 2002; Girolami and
Calderhead,, 2011). We consider MALA updates for the model parameters
$\\{\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}}\\}$ for the mean
model. The dispersion model coefficients are sampled depending on the choice
of likelihood, i.e., $\gamma$ is sampled using a MALA if a saddle-point
approximation of the likelihood is considered, otherwise $\gamma$ is sampled
using a MALA with a numerical approximation to the derivative of the
conditional posterior for $\gamma$ or using a M-H random walk. The parameters
$\\{\xi,\phi_{s}\\}$ are updated using a M-H random walk. All the other
remaining parameters are sampled using Gibbs sampling. In particular, we
employ block updates for
$\mbox{\boldmath$\beta$}_{\mbox{\boldmath{$\mathbf{w}$}}}=\\{\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}}\\}$
and $\gamma$. Proposal variances feature adaptive scaling such that the
optimal acceptance rate ($\approx 58\%$) to capture Langevin dynamics is
achieved upon convergence (see, Carlin and Louis,, 2008; Girolami and
Calderhead,, 2011). Proposal variances in the M-H updates also feature
adaptive scaling such that the optimal acceptance rate ($\approx 33\%$) for
random walks is achieved upon convergence. We outline the full sampling
algorithm at the end of Section 7 of the Supplement. For the hierarchical DGLM
in (4), the specification of a spatial effect translates to fitting a random
intercept mean model. Consequently, having an additional overall intercept
$\beta_{0}$ in the model renders it unidentifiable (see, Gelfand et al.,,
1995, 1996). Hence, $\beta_{0}$ is not estimable, although
$\beta_{0}+\mbox{\boldmath{$\mathbf{w}$}}$ is estimable. $\beta_{0}$ is
estimated through hierarchical centering of the posterior for $\mathbf{w}$
(see, for e.g., Gelfand et al.,, 1996).

The MCMC samples of $\beta$ and $\gamma$ explore their conditional posterior
distributions and point estimates for these model parameters can be obtained
using maximum a-posteriori (MAP) estimates or the posterior means. Although we
obtain point estimates, these estimates do not yield exact zero values since
we have considered a continuous spike-and-slab prior with a spike at $\nu_{0}$
(a small positive number). Additionally, these estimates do not make use of
the all the MCMC samples. We use a Bayesian model averaging–based strategy
that leverages all the MCMC samples to build inference (Hoeting et al.,,
1999). Specifically, we use a FDR-based using Bayesian model averaging
combined with point estimates (see for e.g., Morris et al.,, 2008; Mohammed et
al.,, 2021). Let $\beta_{u}^{(m)}$ for $m=1,\ldots,M$ denote the MCMC samples
(after burn-in and thinning) for the coefficients of the mean model. We
compute $p_{u}=\frac{1}{M}\sum_{m}I\big{(}|\beta_{u}^{(m)}|\leq c\big{)}$,
where $I(\cdot)$ is the indicator function; these probabilities $p_{u}$ can be
interpreted as local FDR (Morris et al.,, 2008). The probability $(1-p_{u})$
can be interpreted as the probability that covariate $u$ is significantly
associated with the response. We use $p_{u}$s to decide on which covariates to
select while controlling the FDR at level $\alpha$. That is, we infer that the
covariate $u$ has a non-zero coefficient if $p_{u}<\kappa_{\alpha}$ for some
threshold $\kappa_{\alpha}\in(0,1)$. We compute the threshold
$\kappa_{\alpha}$ as follows: We first sort the probabilities $p_{u}$ and
denote the sorted probabilities as $p_{(u)}$ for $u=1,\ldots,p$. We then
assign $\kappa_{\alpha}=p_{(u^{*})}$, where
$u^{*}=\max\left\\{\widetilde{u}\mid\frac{1}{\widetilde{u}}\sum_{u=1}^{\widetilde{u}}p_{(u)}\leq\alpha\right\\}$.
This approach caps our false discoveries of selected variables at
$100\alpha\%$. We employ the same approach using the MCMC samples of $\gamma$
to select significant coefficients for the dispersion models.

Posterior inference on $\beta$, $\gamma$ is performed using MAP (point)
estimates along with posterior mean, median, standard deviation and highest
posterior density (HPD) intervals. For $\mathbf{w}$, we employ posterior mean,
median, standard deviation and HPD intervals to perform inference. Next, we
demonstrate some synthetic experiments that document the performance of our
proposed models. The computation has been performed in the R statistical
environment. The required subroutines can be accessed via an open-source
repository at: https://github.com/arh926/sptwdglm.

Table 2: Parameter settings used to obtain varying proportion of zeros in the synthetic data. Proportion of 0s | $\mu_{\mbox{\boldmath$\beta$}_{-0}}\left(\sigma_{\mbox{\boldmath$\beta$}_{-0}}\right)$ | $\gamma_{0}$ | $\mu_{\mbox{\boldmath$\gamma$}_{-0}}\left(\sigma_{\mbox{\boldmath$\gamma$}_{-0}}\right)$
---|---|---|---
15% | 0.50 (0.1) | -1.50 | 0.50 (0.1)
30% | 0.50 (0.1) | 0.70 | 0.50 (0.1)
60% | 0.50 (0.1) | 2.50 | 0.50 (0.1)
80% | 1.00 (0.1) | 4.50 | 0.50 (0.1)
95% | 1.00 (0.1) | 7.00 | 0.50 (0.1)

##  3 Synthetic Experiments

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 $\theta$ 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
$\beta$’s and $\gamma$’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
$\mathbf{X}$ and $\mathbf{Z}$ are hierarchically centered and scaled,
independently sampled Gaussian variables with mean 0 and variance 1.
Naturally, specification of the true $\beta$, $\mathbf{w}$ and $\gamma$
parameters determine the proportion of zeros in the synthetic response. Table
2 contains the parameter specifications used. The true value of the index
parameter, $\xi=1.5$. In an attempt to produce a synthetic setup that
resembles reality we simulate, $\mbox{\boldmath{$\mathbf{w}$}}\sim
N(5(\sin(3\pi\mbox{\boldmath{$\mathbf{s}$}}_{1})+\cos(3\pi\mbox{\boldmath{$\mathbf{s}$}}_{2})),1)$
(see Figure 1, second row). The alternative route would be to fix values of
$\sigma^{2}$ and $\phi_{s}$ and generate a realization
$\mbox{\boldmath{$\mathbf{w}$}}\sim
N_{L}(\mbox{\boldmath{$\mathbf{0}$}}_{L},\sigma^{2}\mbox{\boldmath{$\mathbf{R}$}}(\phi_{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_{\xi}=1$, $b_{\xi}=2$,
$a_{\sigma}=a_{\sigma_{\beta}}=a_{\sigma_{\gamma}}=2$,
$b_{\sigma}=b_{\sigma_{\beta}}=b_{\sigma_{\gamma}}=1$ (producing inverse-Gamma
priors with mean 1 and infinite variance), $a_{\phi_{s}}=0$,
$b_{\phi_{s}}=30$, $\sigma^{-2}_{\beta}=\sigma^{-2}_{\gamma}=10^{-6}$ and
$\nu_{0}=5\times 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\times 10^{3},5\times
10^{3},1\times 10^{4}\\}$ and the number of locations are $L=1\times 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.
$MSE(\mbox{\boldmath$\beta$})=\frac{1}{p}\sum_{i_{\beta}=1}^{p}(\beta_{i_{\beta}}-\widehat{\beta}_{i_{\beta}})^{2}$,
which can be computed similarly for the other parameters. We also compute
average coverage probabilities, for e.g. considering these probabilities for
$\beta$ we define
$CP(\mbox{\boldmath$\beta$})=\frac{1}{M}\sum_{m=1}^{M}I(\mbox{\boldmath$\beta$}_{true}\in(l_{m}(\mbox{\boldmath$\beta$}),u_{m}(\mbox{\boldmath$\beta$})))$,
where $l_{m}(\mbox{\boldmath$\beta$})$ and $u_{m}(\mbox{\boldmath$\beta$})$
are the lower and upper 95% HPD intervals respectively for $\beta$ in
replication $m$; we obtain coverage probabilities for $\mathbf{w}$ and
$\gamma$ 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.

Figure 1: Plots showing synthetic spatial patterns, pattern 1 (top, left
column) and pattern 2 (bottom, left column) and corresponding logarithm of
aggregated synthetic response (right column). Figure 2: Results for synthetic
experiments showing model performance, where MSE is plotted against proportion
of zeros in the synthetic response which is tabulated across models (M1, M2,
M3, M4) vs.
$\\{\mbox{\boldmath$\theta$}_{m},\mbox{\boldmath$\theta$}_{pr}\\}\times$
configuration. Table 3: Table showing results of synthetic experiments for
model selection for models M2 and M4. Corresponding standard deviations are
shown in brackets below.

$N$ | True Overlap | Prop. of 0s | M2 | M4
---|---|---|---|---
Overlap | FPR | TPR | Overlap | FPR | TPR
5000 | 0.00 | 0.15 | 0.03 | 0.04 | 0.69 | 0.00 | 0.00 | 0.90
(0.06) | (0.06) | (0.07) | (0.00) | (0.00) | (0.00)
0.30 | 0.05 | 0.05 | 0.89 | 0.00 | 0.00 | 1.00
(0.09) | (0.09) | (0.10) | (0.00) | (0.00) | (0.00)
0.60 | 0.01 | 0.01 | 0.98 | 0.00 | 0.00 | 1.00
(0.04) | (0.04) | (0.04) | (0.00) | (0.00) | (0.00)
0.80 | 0.04 | 0.04 | 0.99 | 0.00 | 0.00 | 1.00
(0.06) | (0.06) | (0.03) | (0.00) | (0.00) | (0.00)
| 0.95 | 0.10 | 0.06 | 0.89 | 0.08 | 0.04 | 0.95
| (0.20) | (0.10) | (0.20) | (0.10) | (0.11) | (0.08)
0.50 | 0.15 | 0.24 | 0.01 | 0.67 | 0.50 | 0.00 | 0.90
(0.16) | (0.05) | (0.09) | (0.00) | (0.00) | (0.03)
0.30 | 0.49 | 0.01 | 0.94 | 0.50 | 0.00 | 1.00
(0.08) | (0.05) | (0.06) | (0.00) | (0.00) | (0.00)
0.60 | 0.51 | 0.00 | 0.93 | 0.50 | 0.00 | 1.00
(0.03) | (0.00) | (0.06) | (0.00) | (0.00) | (0.00)
0.80 | 0.56 | 0.07 | 0.97 | 0.49 | 0.01 | 1.00
(0.09) | (0.08) | (0.04) | (0.02) | (0.05) | (0.00)
| 0.95 | 0.65 | 0.20 | 0.85 | 0.48 | 0.04 | 0.92
| (0.31) | (0.17) | (0.22) | (0.10) | (0.15) | (0.18)
1.00 | 0.15 | 0.40 | 0.03 | 0.66 | 1.00 | 0.00 | 0.90
(0.25) | (0.05) | (0.07) | (0.00) | (0.00) | (0.00)
0.30 | 0.90 | 0.04 | 0.86 | 1.00 | 0.00 | 1.00
(0.14) | (0.08) | (0.08) | (0.00) | (0.00) | (0.00)
| 0.60 | 1.00 | 0.00 | 0.96 | 1.00 | 0.00 | 1.00
| (0.00) | (0.00) | (0.05) | (0.00) | (0.00) | (0.00)
| 0.80 | 1.00 | 0.00 | 1.00 | 1.00 | 0.00 | 1.00
| (0.00) | (0.00) | (0.00) | (0.00) | (0.00) | (0.00)
|  | 0.95 | 0.35 | 0.15 | 0.85 | 0.89 | 0.10 | 0.89
|  | (0.21) | (0.16) | (0.22) | (0.21) | (0.11) | (0.10)

From the results shown, we see that models M1 and M2 perform poorly.
Estimates, $\widehat{\mbox{\boldmath$\beta$}}$ remain fairly unaffected as
compared to $\widehat{\mbox{\boldmath$\gamma$}}$ and $\widehat{\xi}$, 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
$\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}}$ and $\gamma$. Chain
lengths are set to $1\times 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, $\widehat{\mbox{\boldmath$\theta$}}$ 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 ($\approx 1$) across all settings; only declining marginally
for $\mathbf{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.

##  4 Supplementary Analysis

The online Supplement to this paper contains details on the derivations of the
posteriors essential for constructing MCMC subroutines. They are outlined in
Section 7. Section 8 features additional simulation experiments that
supplement those outlined previously in Section 3. It documents performance of
M2, shown in Table 7, contains results of experiments for scenarios featuring
spatial covariates, shown in Table 8, and varying spatial patterns as seen in
Figure 1, shown in Tables 9 and 10. Convergence diagnostics are shown for
selected model parameters (index parameter $\xi$) in Section 8.2. Contents of
the R-package are described in Section 8.3. Finally, results for models M1 and
M3 pertaining to the real data analysis described in the next section appear
in Section 9, Tables 11, 12, 13 and 14.

##  5 Automobile Insurance Premiums, Connecticut, 2008

Figure 3: (left) Spatial plot of zip-code level aggregated pure-premium
$\times 10^{-6}$ for the state of Connecticut, 2008 (right) histogram for the
pure-premium overlaid with a probability density estimate.

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.

  * b.

Policy level: (i) policy payments, measured in United States dollars, (ii)
exposure–measured in policy years, for e.g. 0.5 indicates a coverage period of
6 months or, half a year, measured in years, (iii) policy risk–having two
levels, which is assigned by the insurance company based on provided
information by the individual, (iv) deductible limit–with 8 categories.

  * c.

Spatial: 5-digit zip code.

Derived variables like age categories, vehicle age in years and interactions
like gender $\times$ marital status are computed and used as covariates in the
model. For the state of Connecticut, 1,513,655 ($\approx$ 1.5 million) data
records were obtained in the year 2008, at 281 zip-codes. Zip-codes are areal
units, we consider the latitude-longitude corresponding to the centroid of
each zip code as the point reference counterpart unit for our application
purposes. Distance between two zip-codes is then specified as the Euclidean
distance between their centroids. The proportion of zeros in the payments is
95.73%. From an insurer’s perspective, policy rate-making is the problem of
assigning policy-premium to a new customer’s policy based on their covariate
information (for instance, individual level, policy and residence zip-code).
We achieve this via out-sample prediction. To that end, we consider a 60-40
split for the data, the split is performed by using stratified sampling
without replacement over zip-codes, such that the same 281 zip-codes are also
available in the out-sample. The training data then contains $N_{tr}=908,741$
observations with $N_{pr}=604,914$ observations kept in reserve for prediction
constituting the out-sample data.

We denote payments towards a policy made by individual $i$, residing in zip-
code $j$ as $y_{ij}$ with an exposure of $t_{ij}$. We assume that the policy-
premium defined as $y^{*}_{ij}=\frac{y_{ij}}{t_{ij}}\sim
Tw\left(\mu_{ij},\phi_{ij},\xi\right)$, which implies $y_{ij}\sim
Tw(t_{ij}\mu_{ij},t_{ij}^{2-\xi}\phi_{ij},\xi)$ using the scale invariance
property. The following hierarchical DGLM is then specified,

$\displaystyle\begin{split}\log\mu_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})&=-\log
t_{ij}+\mbox{\boldmath{$\mathbf{x}$}}_{ij}^{T}(\mbox{\boldmath{$\mathbf{s}$}}_{i})\mbox{\boldmath$\beta$}+\mbox{\boldmath{$\mathbf{f}$}}_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})^{T}\mbox{\boldmath{$\mathbf{w}$}}(\mbox{\boldmath{$\mathbf{s}$}}_{i}),\\\
\log\phi_{ij}&=-(2-\xi)\log
t_{ij}+\mbox{\boldmath{$\mathbf{z}$}}_{ij}^{T}\mbox{\boldmath$\gamma$},\end{split}$
(7)

when considered with a spatial specification, where the terms $-\log t_{ij}$
and $-(2-\xi)\log t_{ij}$ act as offsets for the respective mean and
dispersion models. Given the covariates described at the beginning $p=q=29$,
producing a $N_{tr}\times(p-1)$ design matrix for the mean model and a
$N_{tr}\times q$ design matrix for the dispersion model. The model in (7)
specifies model M3 from Table 1, M4 is obtained by specifying
$\pi(\mbox{\boldmath$\theta$}_{vs})$ from (6) on
$\mbox{\boldmath$\theta$}_{m}$. M1 is obtained by setting
$f_{ij}(\mbox{\boldmath{$\mathbf{s}$}}_{i})=\mbox{\boldmath{$\mathbf{0}$}}$
and M2 is obtained by specifying $\pi(\mbox{\boldmath$\theta$}_{vs})$ on
$\mbox{\boldmath$\theta$}_{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 from fitting models M2 and M4, controlling for FDR at
1%. The performance of M1–4 is assessed using the Akaike Information Criteria
(AIC) (Akaike,, 1998).

Figure 4: Spatial plot showing the posterior estimates for 281 zipcodes from
(left) M3, and (right) M4 in Connecticut. The locations are color coded based
on significance, with white indicating a location with 0 in its HPD interval,
blue (red) indicating HPD interval with both endpoints negative (positive).

While specifying (5) and (6) the hyper-parameter settings used are,
$\sigma^{2}_{\beta}=\sigma^{2}_{\gamma}=10^{6}$,
$a_{\sigma_{\beta}}=a_{\sigma}=a_{\sigma_{\gamma}}=2$ and
$b_{\sigma_{\beta}}=b_{\sigma}=b_{\sigma_{\gamma}}=1$ (generating inverse-
gamma priors with mean 1 and infinite variance), $a_{\phi_{s}}=0$,
$b_{\phi_{s}}=60$. We maintain $a_{\xi}=1$ and $b_{\xi}=2$ for all models. For
ease of implementation, the fractal parameter, $\nu$ is fixed at 0.5,
producing the exponential covariance kernel. We consider $1\times 10^{5}$ MCMC
iterations for generating samples from respective posteriors with burn-in
diagnosed at $5\times 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 11, 12, 13 and 14 in the Supplementary
Materials. Posterior estimates for the spatial effects in models M3 and M4 are
shown in Figure 4. Zip-codes with significant effects are color coded
appropriately.

Table 4: Estimated Coefficients for fixed effects corresponding to model M2.
We show the maximum a-posteriori (MAP) estimates along with median, mean,
standard deviation and highest posterior density intervals (HPD). The
variables listed are a result of FDR-based variable selection at 1%.

Parameters | Levels | MAP | Median | Mean | SD | Lower HPD | Upper HPD
---|---|---|---|---|---|---|---
$\beta$ | (Intercept) | – | 5.815 | 5.965 | 5.955 | 0.151 | 5.706 | 6.227
age.car | 1 | 0.113 | 0.112 | 0.110 | 0.033 | 0.044 | 0.173
2 | 0.234 | 0.239 | 0.240 | 0.030 | 0.186 | 0.304
5 | 0.124 | 0.122 | 0.121 | 0.033 | 0.062 | 0.189
risk | S | -0.213 | -0.217 | -0.217 | 0.023 | -0.258 | -0.168
agec | 3 | -0.399 | -0.399 | -0.400 | 0.030 | -0.454 | -0.341
4 | -0.518 | -0.515 | -0.515 | 0.035 | -0.580 | -0.444
5 | -0.712 | -0.710 | -0.708 | 0.065 | -0.832 | -0.575
gender | F | 0.493 | 0.502 | 0.513 | 0.067 | 0.400 | 0.648
M | 0.608 | 0.614 | 0.619 | 0.053 | 0.521 | 0.731
marital | M | -0.376 | -0.392 | -0.403 | 0.067 | -0.577 | -0.302
deductible | B | 0.898 | 1.051 | 1.104 | 0.244 | 0.740 | 1.525
E | 0.616 | 0.482 | 0.478 | 0.158 | 0.192 | 0.741
F | 0.580 | 0.428 | 0.435 | 0.154 | 0.135 | 0.656
| G | 0.287 | 0.348 | 0.358 | 0.159 | 0.073 | 0.633
$\gamma$ | (Intercept) | – | 7.354 | 7.346 | 7.345 | 0.048 | 7.249 | 7.435
age.car | 0 | -0.820 | -0.811 | -0.811 | 0.053 | -0.911 | -0.714
1 | -1.024 | -1.017 | -1.016 | 0.051 | -1.115 | -0.922
2 | -0.888 | -0.874 | -0.873 | 0.052 | -0.969 | -0.776
3 | -0.864 | -0.854 | -0.851 | 0.052 | -0.953 | -0.760
4 | -0.843 | -0.838 | -0.838 | 0.052 | -0.936 | -0.743
5 | -0.788 | -0.781 | -0.780 | 0.052 | -0.877 | -0.682
6 | -0.770 | -0.765 | -0.765 | 0.053 | -0.862 | -0.665
7 | -0.730 | -0.723 | -0.723 | 0.053 | -0.829 | -0.627
risk | S | 0.114 | 0.114 | 0.114 | 0.011 | 0.093 | 0.135
agec | 5 | -0.289 | -0.293 | -0.293 | 0.029 | -0.345 | -0.234
deductible | B | -0.971 | -0.963 | -0.951 | 0.099 | -1.136 | -0.759
C | -0.519 | -0.522 | -0.517 | 0.049 | -0.614 | -0.414
D | -0.565 | -0.568 | -0.564 | 0.044 | -0.657 | -0.470
E | -0.506 | -0.501 | -0.495 | 0.042 | -0.568 | -0.394
F | -0.348 | -0.346 | -0.340 | 0.039 | -0.410 | -0.242
H | 0.411 | 0.419 | 0.426 | 0.097 | 0.251 | 0.618
genderMarital | A | -0.113 | -0.106 | -0.097 | 0.043 | -0.168 | -0.001
$\xi$ | – | – | 1.673 | 1.673 | 1.673 | 0.001 | 1.671 | 1.676

Table 5: Estimated coefficients for fixed effects corresponding to model M4.
We show the MAP, median, mean, standard deviation and HPDs. The variables
listed are a result of FDR-based variable selection at 1%.

Parameters | Levels | MAP | Median | Mean | SD | Lower HPD | Upper HPD
---|---|---|---|---|---|---|---
$\beta$ | (Intercept) | – | 5.215 | 5.213 | 5.216 | 0.135 | 5.165 | 5.327
age.car | 0 | 0.252 | 0.251 | 0.251 | 0.063 | 0.123 | 0.377
1 | 0.329 | 0.317 | 0.318 | 0.058 | 0.217 | 0.431
2 | 0.424 | 0.423 | 0.425 | 0.061 | 0.318 | 0.549
3 | 0.209 | 0.193 | 0.195 | 0.060 | 0.084 | 0.318
4 | 0.236 | 0.252 | 0.255 | 0.062 | 0.149 | 0.387
5 | 0.297 | 0.302 | 0.307 | 0.059 | 0.191 | 0.428
6 | 0.182 | 0.209 | 0.214 | 0.061 | 0.113 | 0.343
7 | 0.155 | 0.154 | 0.157 | 0.063 | 0.049 | 0.286
risk | S | -0.190 | -0.184 | -0.183 | 0.023 | -0.226 | -0.137
agec | 2 | -0.121 | -0.126 | -0.128 | 0.030 | -0.186 | -0.068
3 | -0.461 | -0.467 | -0.470 | 0.031 | -0.531 | -0.412
4 | -0.608 | -0.617 | -0.619 | 0.034 | -0.685 | -0.554
5 | -0.721 | -0.699 | -0.696 | 0.065 | -0.808 | -0.561
gender | F | 0.535 | 0.540 | 0.541 | 0.066 | 0.407 | 0.679
M | 0.670 | 0.705 | 0.720 | 0.100 | 0.535 | 0.919
deductible | B | 1.337 | 1.383 | 1.408 | 0.216 | 1.036 | 1.866
C | 0.209 | 0.256 | 0.303 | 0.152 | 0.073 | 0.624
D | 0.603 | 0.627 | 0.672 | 0.161 | 0.397 | 0.982
E | 0.800 | 0.825 | 0.864 | 0.157 | 0.625 | 1.189
F | 0.753 | 0.774 | 0.819 | 0.155 | 0.592 | 1.154
G | 0.756 | 0.786 | 0.825 | 0.158 | 0.580 | 1.169
H | 0.820 | 0.829 | 0.824 | 0.190 | 0.463 | 1.149
| genderMarital | B | -0.368 | -0.236 | -0.169 | 0.233 | -0.457 | 0.311
$\gamma$ | (Intercept) | – | 6.423 | 6.429 | 6.415 | 0.073 | 6.310 | 6.504
age.car | 0 | -0.790 | -0.809 | -0.811 | 0.040 | -0.891 | -0.731
1 | -1.006 | -1.025 | -1.028 | 0.039 | -1.099 | -0.954
2 | -0.875 | -0.885 | -0.889 | 0.039 | -0.972 | -0.819
3 | -0.842 | -0.861 | -0.863 | 0.039 | -0.950 | -0.799
4 | -0.831 | -0.844 | -0.848 | 0.039 | -0.923 | -0.777
5 | -0.781 | -0.792 | -0.795 | 0.039 | -0.879 | -0.728
6 | -0.762 | -0.766 | -0.769 | 0.039 | -0.849 | -0.703
7 | -0.725 | -0.732 | -0.735 | 0.039 | -0.814 | -0.662
risk | S | 0.110 | 0.110 | 0.110 | 0.012 | 0.087 | 0.132
agec | 5 | -0.262 | -0.262 | -0.261 | 0.030 | -0.316 | -0.201
deductible | B | -1.023 | -1.026 | -1.029 | 0.158 | -1.328 | -0.721
C | -0.534 | -0.553 | -0.571 | 0.069 | -0.700 | -0.461
D | -0.591 | -0.611 | -0.634 | 0.067 | -0.760 | -0.534
E | -0.533 | -0.552 | -0.575 | 0.066 | -0.702 | -0.478
F | -0.377 | -0.390 | -0.416 | 0.065 | -0.543 | -0.335
H | 0.352 | 0.361 | 0.369 | 0.107 | 0.170 | 0.591
$\xi$ | – |  | 1.667 | 1.667 | 1.667 | 0.001 | 1.665 | 1.670

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, Halder et al.,, 2021) but within a more robust
proposed model choice framework. This is evident when comparing model
estimates between Table 5 and Table 13 and 14 in the Supplement. The marital
status and interaction between gender and marital status not being selected.
We conclude by observing that the spatial effects are significantly positive
in major cities in Connecticut indicating higher spatial risk as opposed to
sparsely populated regions showing significantly lower risk.

Table 6: Table showing AIC and out-sample square root deviance for models
M1–M4.

| M1 | M2 | M3 | M4
---|---|---|---|---
AIC | 1340060 | 1117733 | 1117114 | 1115363
$\sqrt{\mbox{ Deviance }}$ | 5565.209 | 5509.549 | 5507.926 | 5441.23

##  6 Discussion

Double generalized linear models have not seen much use after their inception
by Lee and Nelder, (2006). 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., Heaton et al.,, 2019) 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, Cressie,, 2015, 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.

Supplementary Materials for
“Bayesian Variable Selection in Double Generalized Linear Tweedie Spatial
Process Models”

##  7 Posteriors

Under chosen priors on
$\\{\mbox{\boldmath$\theta$}_{m},\mbox{\boldmath$\theta$}_{pr}\\}=\\{\mbox{\boldmath$\beta$},\mbox{\boldmath$\gamma$},\xi,\sigma^{2},\phi_{s}\\}$,
the resulting joint posterior is specified by

$\displaystyle\begin{aligned}
\pi\Big{(}\mbox{\boldmath$\theta$}_{m},\mbox{\boldmath$\theta$}_{pr}\mid\mbox{\boldmath{$\mathbf{y}$}}\Big{)}&\propto
U\Big{(}\xi\mid a_{\xi},b_{\xi}\Big{)}\times
N_{p}\Big{(}\mbox{\boldmath$\beta$}\mid\mbox{\boldmath{$\mathbf{0}$}}_{p},\mbox{\boldmath$\lambda$}_{\beta}^{T}\mbox{\boldmath{$\mathbf{I}$}}_{p}\Big{)}\times
N_{q}\Big{(}\mbox{\boldmath$\gamma$}\mid\mbox{\boldmath{$\mathbf{0}$}}_{q},\mbox{\boldmath$\lambda$}_{\gamma}^{T}\mbox{\boldmath{$\mathbf{I}$}}_{q}\Big{)}\times\\\
&\hskip 14.22636ptU\Big{(}\phi_{s}\mid
a_{\phi_{s}},b_{\phi_{s}}\Big{)}\times{\rm Gamma}\Big{(}\sigma^{-2}\mid
a_{\sigma},b_{\sigma}\Big{)}\times
N_{L}\Big{(}\mbox{\boldmath{$\mathbf{w}$}}\mid\mbox{\boldmath{$\mathbf{0}$}}_{L},\sigma^{2}\mbox{\boldmath{$\mathbf{R}$}}(\phi_{s})\Big{)}\times\\\
&\hskip
14.22636pt{Tw}\Big{(}\mbox{\boldmath{$\mathbf{y}$}}\mid\mbox{\boldmath$\mu$}=\exp(\mbox{\boldmath{$\mathbf{X}$}}\mbox{\boldmath$\beta$}+\mbox{\boldmath{$\mathbf{F}$}}\mbox{\boldmath{$\mathbf{w}$}}),\mbox{\boldmath$\phi$}=\exp(\mbox{\boldmath{$\mathbf{Z}$}}\mbox{\boldmath$\gamma$}),\xi\Big{)},\end{aligned}$
(8)

We list the posteriors for individual parameters in the following equations,

$\displaystyle\pi\left(\mbox{\boldmath$\beta$}\mid\mbox{\boldmath$\gamma$},\xi\right)\propto\exp\left\\{-\left(\sum\limits_{i=1}^{L}\sum\limits_{j=1}^{n_{i}}\frac{1}{\phi_{ij}}d(y_{ij}\mid\mu_{ij}(\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}}),\xi)+\frac{\sigma^{-2}_{\beta}}{2}||\mbox{\boldmath$\beta$}||_{2}^{2}+\frac{\sigma^{-2}}{2}\mbox{\boldmath{$\mathbf{w}$}}^{T}\mbox{\boldmath{$\mathbf{R}$}}^{-1}(\phi_{s})\mbox{\boldmath{$\mathbf{w}$}}\right)\right\\},$
(9)
$\displaystyle\pi\left(\mbox{\boldmath$\gamma$}\mid\mbox{\boldmath$\beta$},\xi\right)\propto\exp\left\\{-\left(\sum\limits_{i=1}^{L}\sum\limits_{j=1}^{n_{i}}\frac{1}{\phi_{ij}(\mbox{\boldmath$\gamma$})}d(y_{ij}\mid\mu_{ij},\xi)+\frac{\log\phi_{ij}(\mbox{\boldmath$\gamma$})}{2}I(y_{ij}>0)+\frac{\sigma^{-2}_{\gamma}}{2}||\mbox{\boldmath$\gamma$}||_{2}^{2}\right)\right\\},$
$\displaystyle\pi(\xi\mid\mbox{\boldmath$\beta$},\mbox{\boldmath$\gamma$})\propto\prod_{i=1}^{L}\prod_{j=1}^{n_{i}}c_{ij}(y_{ij}\mid\phi_{ij},\xi)\exp\left\\{-\frac{1}{\phi_{ij}}d(y_{ij}\mid\mu_{ij},\xi)\right\\}I\left(\xi\in(a_{\xi},b_{\xi})\right),$
$\displaystyle\pi(\phi_{s})\propto|\mbox{\boldmath{$\mathbf{R}$}}(\phi_{s})|^{-1/2}\exp\left(-\frac{\sigma^{-2}}{2}\mbox{\boldmath{$\mathbf{w}$}}^{T}\mbox{\boldmath{$\mathbf{R}$}}^{-1}(\phi_{s})\mbox{\boldmath{$\mathbf{w}$}}\right)I\left(\phi_{s}\in(a_{\phi_{s}},b_{\phi_{s}})\right),$
$\displaystyle\pi(\sigma^{-2})={\rm
Gamma}\left(a_{\sigma}+\frac{L}{2},b_{\sigma}+\frac{1}{2}\mbox{\boldmath{$\mathbf{w}$}}^{T}\mbox{\boldmath{$\mathbf{R}$}}^{-1}(\phi_{s})\mbox{\boldmath{$\mathbf{w}$}}\right),$

where $|\cdot|$ denotes the determinant and $d(y_{ij}\mid\mu_{ij},\xi)$ is the
deviance function. In the scenario where we have no spatial effect in the
DGLM, the reduced set of parameters are
$\mbox{\boldmath$\theta$}_{m}=\\{\mbox{\boldmath$\beta$},\mbox{\boldmath$\gamma$},\xi\\}$
having a similar joint posterior after omitting the second row involving
spatial process prior specification and setting
$\mbox{\boldmath$\mu$}=\exp(\mbox{\boldmath{$\mathbf{X}$}}\mbox{\boldmath$\beta$})$
in the likelihood, keeping the same prior specifications on the other
parameters, the resulting posteriors are as follows,

$\displaystyle\pi\left(\mbox{\boldmath$\beta$}\mid\mbox{\boldmath$\gamma$},\xi\right)\propto\exp\left\\{-\left(\sum\limits_{i=1}^{N}\frac{1}{\phi_{i}}d(y_{i}\mid\mu_{i}(\mbox{\boldmath$\beta$}),\xi)+\frac{\sigma^{-2}_{\beta}}{2}||\mbox{\boldmath$\beta$}||_{2}^{2}\right)\right\\},$
$\displaystyle\pi\left(\mbox{\boldmath$\gamma$}\mid\mbox{\boldmath$\beta$},\xi\right)\propto\exp\left\\{-\left(\sum\limits_{i=1}^{N}\frac{1}{\phi_{i}(\mbox{\boldmath$\gamma$})}d(y_{i}\mid\mu_{i},\xi)+\frac{1}{2}\log\phi_{ij}(\mbox{\boldmath$\gamma$})I(y_{ij}>0)+\frac{\sigma^{-2}_{\gamma}}{2}||\mbox{\boldmath$\gamma$}||_{2}^{2}\right)\right\\},$
$\displaystyle\pi(\xi)\propto\prod_{i=1}^{N}c_{i}(y_{i}\mid\phi_{i},\xi)\exp\left\\{-\frac{1}{\phi_{i}}d(y_{i}\mid\mu_{i},\xi)\right\\}I\left(\xi\in(a_{\xi},b_{\xi})\right).$

Updates leveraging MALA for
${\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}}}$ (or $\beta$) and
$\gamma$ require the proposals to be specified appropriately using the
gradients of log-posterior densities,
$\nabla\log\pi\left(\mbox{\boldmath$\beta$}\mid\mbox{\boldmath$\gamma$},\xi\right)$
and
$\nabla\log\pi\left(\mbox{\boldmath$\gamma$}\mid\mbox{\boldmath$\beta$},\xi\right)$
respectively. Candidate samples are obtained using,

$\displaystyle{(\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}})^{T}}^{*}=(\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}})^{T}+\frac{\tau_{\beta,w}^{2}}{2}\mbox{\boldmath{$\mathbf{A}$}}_{\beta,w}\nabla\log\pi\left(\mbox{\boldmath$\beta$}\mid\mbox{\boldmath$\gamma$},\xi\right)+\tau_{\beta,w}\mbox{\boldmath{$\mathbf{A}$}}_{\beta,w}^{1/2}\cdot
N_{p+L}(\mbox{\boldmath{$\mathbf{0}$}},\mbox{\boldmath{$\mathbf{I}$}}_{p+L}),$
(10)
$\displaystyle\mbox{\boldmath$\gamma$}^{*}=\mbox{\boldmath$\gamma$}+\frac{\tau_{\gamma}^{2}}{2}\mbox{\boldmath{$\mathbf{A}$}}_{\gamma}\nabla\log\pi\left(\gamma\mid\mbox{\boldmath$\beta$},\xi\right)+\tau_{\gamma}\mbox{\boldmath{$\mathbf{A}$}}_{\gamma}^{1/2}\cdot
N_{q}(\mbox{\boldmath{$\mathbf{0}$}},\mbox{\boldmath{$\mathbf{I}$}}_{q}),$

where $\mbox{\boldmath{$\mathbf{A}$}}_{\beta,w}^{-1}={\rm
E}\left[-\nabla^{2}\log\pi\left(\mbox{\boldmath$\beta$}_{\mbox{\boldmath{$\mathbf{w}$}}}\mid-\right)\right]$
and $\mbox{\boldmath{$\mathbf{A}$}}_{\gamma}^{-1}={\rm
E}\left[-\nabla^{2}\log\pi\left(\mbox{\boldmath$\gamma$}\mid-\right)\right]$.

Under a continuous spike and slab prior on $\beta$ and $\gamma$ the posteriors
for the hierarchical latent DGLM is as follows,

$\displaystyle\pi\left(\mbox{\boldmath$\beta$}_{\mbox{\boldmath{$\mathbf{w}$}}}\mid-\right)\propto\exp\left\\{-\left(\sum\limits_{i=1}^{L}\sum\limits_{j=1}^{n_{i}}\frac{1}{\phi_{ij}}d(y_{ij}\mid\mu_{ij}(\mbox{\boldmath$\beta$},\mbox{\boldmath{$\mathbf{w}$}}),\xi)+\frac{1}{2}\mbox{\boldmath$\beta$}^{T}\mbox{\boldmath$\Gamma$}_{\beta}^{-1}\mbox{\boldmath$\beta$}+\frac{\sigma^{-2}}{2}\mbox{\boldmath{$\mathbf{w}$}}^{T}\mbox{\boldmath{$\mathbf{R}$}}^{-1}(\phi_{s})\mbox{\boldmath{$\mathbf{w}$}}\right)\right\\}$
(11)
$\displaystyle\left.\begin{array}[]{l}\mbox{\boldmath$\Gamma$}_{\beta}={\rm
diag}\\{\zeta_{\beta,i_{b}}\sigma_{\beta,i_{b}}^{2}\\}\\\
\pi(\zeta_{\beta,i_{b}}\mid-)=\frac{\alpha_{1,\beta}}{\alpha_{1,\beta}+\alpha_{2,\beta}}\delta_{\nu_{0}}(\cdot)+\frac{\alpha_{2,\beta}}{\alpha_{1,\beta}+\alpha_{2,\beta}}\delta_{\nu_{1}}(\cdot)\\\
\pi(\sigma^{-2}_{\beta,i_{b}}\mid-)={\rm
Gamma}\left(a_{\sigma}+\frac{1}{2},b_{\sigma}+\frac{\beta^{2}_{i_{b}}}{2\zeta_{\beta,i_{b}}}\right)\\\
\pi(\alpha_{\beta}\mid-)={\rm
Beta}\left(1+\\#\\{i_{b}:\zeta_{\beta,i_{b}}=1\\},1+\\#\\{i_{b}:\zeta_{\beta,i_{b}}=\nu_{0}\\}\right)\end{array}\right\\},i_{b}=1,2,\ldots,p$
$\displaystyle\pi\left(\mbox{\boldmath$\gamma$}\mid-\right)\propto\exp\left\\{-\left(\sum\limits_{i=1}^{L}\sum\limits_{j=1}^{n_{i}}\frac{1}{\phi_{ij}(\mbox{\boldmath$\gamma$})}d(y_{ij}\mid\mu_{ij},\xi)+\frac{1}{2}\log\phi_{ij}(\mbox{\boldmath$\gamma$})I(y_{ij}>0)+\frac{1}{2}\mbox{\boldmath$\gamma$}^{T}\mbox{\boldmath$\Gamma$}_{\gamma}^{-1}\mbox{\boldmath$\gamma$}\right)\right\\},$
$\displaystyle\left.\begin{array}[]{l}\mbox{\boldmath$\Gamma$}_{\gamma}={\rm
diag}\\{\zeta_{\gamma,i_{g}}\sigma_{\gamma,i_{g}}^{2}\\}\\\
\pi(\zeta_{\gamma,i_{g}}\mid-)=\frac{\alpha_{1,\gamma}}{\alpha_{1,\gamma}+\alpha_{2,\gamma}}\delta_{\nu_{0}}(\cdot)+\frac{\alpha_{2,\gamma}}{\alpha_{1,\gamma}+\alpha_{2,\gamma}}\delta_{\nu_{1}}(\cdot)\\\
\pi(\sigma^{-2}_{\gamma,i_{g}}\mid-)={\rm
Gamma}\left(a_{\sigma}+\frac{1}{2},b_{\sigma}+\frac{\gamma^{2}_{i_{g}}}{2\zeta_{\gamma,i_{g}}}\right)\\\
\pi(\alpha_{\gamma}\mid-)={\rm
Beta}\left(1+\\#\\{i_{g}:\zeta_{\gamma,i_{g}}=1\\},1+\\#\\{i_{g}:\zeta_{\gamma,i_{g}}=\nu_{0}\\}\right)\end{array}\right\\},i_{g}=1,2,\ldots,q$
$\displaystyle\pi(\xi\mid-)\propto\prod_{i=1}^{L}\prod_{j=1}^{n_{i}}c_{ij}(y_{ij}\mid\phi_{ij},\xi)\exp\left\\{-\frac{1}{\phi_{ij}}d(y_{ij}\mid\mu_{ij},\xi)\right\\}I\left(\xi\in(a_{\xi},b_{\xi})\right),$
$\displaystyle\pi(\phi_{s}\mid-)\propto|\mbox{\boldmath{$\mathbf{R}$}}(\phi_{s})|^{-1/2}\exp\left(-\frac{\sigma^{-2}}{2}\mbox{\boldmath{$\mathbf{w}$}}^{T}\mbox{\boldmath{$\mathbf{R}$}}^{-1}(\phi_{s})\mbox{\boldmath{$\mathbf{w}$}}\right)I\left(\phi_{s}\in(a_{\phi_{s}},b_{\phi_{s}})\right),$
$\displaystyle\pi(\sigma^{-2}\mid-)={\rm
Gamma}\left(a_{\sigma}+\frac{L}{2},b_{\sigma}+\frac{1}{2}\mbox{\boldmath{$\mathbf{w}$}}^{T}\mbox{\boldmath{$\mathbf{R}$}}^{-1}(\phi_{s})\mbox{\boldmath{$\mathbf{w}$}}\right),$

where,

$\displaystyle\alpha_{1,\beta}=(1-\alpha_{\beta})\nu_{0}^{-1/2}\exp\left\\{-\frac{\beta_{i_{b}}^{2}}{2\nu_{0}\sigma_{\beta,i_{b}}^{2}}\right\\},~{}\alpha_{2,\beta}=\alpha_{\beta}\exp\left\\{-\frac{\beta_{i_{b}}^{2}}{2\sigma^{2}_{\beta,i_{b}}}\right\\},i_{b}=1,2,\ldots,p;$
$\displaystyle\alpha_{1,\gamma}=(1-\alpha_{\gamma})\nu_{0}^{-1/2}\exp\left\\{-\frac{\gamma_{i_{g}}^{2}}{2\nu_{0}\sigma_{\gamma,i_{g}}^{2}}\right\\},~{}\alpha_{2,\gamma}=\alpha_{\gamma}\exp\left\\{-\frac{\gamma_{i_{g}}^{2}}{2\sigma^{2}_{\gamma,i_{g}}}\right\\},i_{g}=1,2,\ldots,q.$

Under no latent specification the resulting posteriors can be obtained
similarly, by omitting updates for the latent process and associated process
parameters/hyper-parameters. The joint posterior is sampled by leveraging the
posteriors in eq. (LABEL:eq:posterior-1b) and employing a similar hybrid
sampling strategy as earlier, with additional Gibbs updates for the spike and
slab parameters,
$\\{\mbox{\boldmath$\zeta$}_{\beta},\mbox{\boldmath$\sigma$}^{2}_{\beta},\alpha_{\beta},\mbox{\boldmath$\zeta$}_{\gamma},\mbox{\boldmath$\sigma$}^{2}_{\gamma},\alpha_{\gamma}\\}$.

##  8 Experiments, Diagnostics and Software

We outline supporting simulation experiments, convergence diagnostics and
details on available software to supplement Section 3 in the manuscript. We
begin with results from additional experiments.

###  8.1 Supporting Experiments

We performed simulation experiments to assess performance of model M2—Table 2
shows the results featuring the same settings as outlined in Section 3, only
without the spatial effects.

Configuration | Prop. of 0s | $CP(\mbox{\boldmath$\theta$})$ | False Positive Rate | True Positive Rate
---|---|---|---|---
Configuration 1 | 0.15 | 1.00 (0.00) | 0.00 (0.00) | 1.00 (0.00)
0.30 | 1.00 (0.00) | 0.00 (0.00) | 1.00 (0.00)
0.60 | 1.00 (0.00) | 0.00 (0.00) | 1.00 (0.00)
0.80 | 1.00 (0.10) | 0.04 (0.06) | 0.99 (0.03)
Configuration 2 | 0.15 | 1.00 (0.00) | 0.01 (0.05) | 0.99 (0.03)
0.30 | 1.00 (0.00) | 0.00 (0.00) | 1.00 (0.00)
0.60 | 1.00 (0.10) | 0.00 (0.00) | 0.98 (0.04)
0.80 | 0.90 (0.10) | 0.01 (0.05) | 0.99 (0.03)
Configuration 3 | 0.15 | 0.90 (0.10) | 0.00 (0.00) | 1.00 (0.00)
0.30 | 1.00 (0.00) | 0.01 (0.04) | 1.00 (0.00)
0.60 | 1.00 (0.00) | 0.00 (0.00) | 0.99 (0.03)
0.80 | 1.00 (0.00) | 0.05 (0.09) | 0.96 (0.05)

Table 7: Results for synthetic experiments corresponding to Table 3 for model
M2 showing average coverage probabilities for the estimated model
coefficients.

In what follows, we outline simulation experiments for the scenario where we
have spatial covariates in our data. Finally, we focus on showing simulations
results for different synthetic spatial patterns.

####  8.1.1 Spatial Covariates

Setting up experiments with spatial covariates we are mindful of possible
endogeneity issues arising from included spatial covariates being correlated
with the true spatial pattern (see, for e.g., Fan and Liao,, 2014). We use the
following settings and true values: $N=1\times 10^{4}$, $L=1\times 10^{2}$,
$\mbox{\boldmath$\beta$}=(1.0,1.5,1\times 10^{-5},1.4,1.1,1\times
10^{-5},2.5)^{T}$, $\mbox{\boldmath$\gamma$}=(1.0,1\times
10^{-5},1.5,1.1,1\times 10^{-5},-2.5,1\times 10^{-5})^{T}$. In $\mathbf{X}$
and $\mathbf{Z}$ we include two spatial covariates in the last two columns:
$\mbox{\boldmath{$\mathbf{x}$}}_{6}=\mbox{\boldmath{$\mathbf{z}$}}_{6}\sim
N(5(\cos(3\pi s_{x})+\sin(3\pi s_{y})),1)$ and
$\mbox{\boldmath{$\mathbf{x}$}}_{7}=\mbox{\boldmath{$\mathbf{z}$}}_{7}\sim
N(2(\cos(\pi s_{x})+\sin(\pi s_{y})),1)$, while the rest were sampled from
standard Gaussian distributions resulting in an average of 50% zeros in the
data. Under this setup, we observe that for the mean model
$\mbox{\boldmath{$\mathbf{x}$}}_{6}$ is significant and
$\mbox{\boldmath{$\mathbf{x}$}}_{7}$ is not while for the dispersion model it
is the other way around. The true spatial effect is simulated using
$\sigma^{2}=1.5$ and $\phi_{s}=3$ from an exponential kernel. The true index
parameter, $\xi=1.5$. The design matrices were centered and scaled. The
resulting absolute value of the correlations within covariates and between
them and the true spatial effect from the setup above was $<0.01$. To assess
the performance of our models we focus on M3 and and M4 and compute the MSE,
coverage probability, FPR and TPR for each of these models. We performed 100
replications under these settings. The results are shown below in Table 8.
Models M3 and M4 were able to estimate coefficients $\beta_{6}$, $\beta_{7}$
and $\gamma_{6}$, $\gamma_{7}$ with sufficient accuracy.

Model | MSE ($\theta$) | MSE ($\mathbf{w}$) | CP ($\theta$) | CP ($\mathbf{w}$) | FPR | TPR
---|---|---|---|---|---|---
M3 | 0.008 | 0.072 | 0.907 | 0.940 | – | –
(0.019) | (0.186) | (0.111) | (0.092) | – | –
M4 | 0.003 | 0.062 | 0.916 | 0.959 | 0.000 | 1.000
(0.004) | (0.129) | (0.063) | (0.09) | (0.165) | (0.035)
Table 8: Table showing results for simulation carried using spatial covariates
in design matrices $\mathbf{X}$ and $\mathbf{Z}$.

####  8.1.2 Spatial Patterns

We show the results for simulation experiments performed for patters shown in
Figure 1 of the manuscript in Tables 9 and 10 below.

$N$ | Prop. of 0s | MSE | CP($\theta$)
---|---|---|---
$\beta$ | $\gamma$ | $\mathbf{w}$ | $\xi$ | $\sigma^{2}$ | $\phi_{s}$ | $\mbox{\boldmath$\theta$}_{m}$ | $\mbox{\boldmath$\theta$}_{pr}$
2000 | 0.15 | 0.00 | 0.00 | 0.02 | 0.00 | 0.04 | 1.45 | 1.00 | 0.95
(0.06) | (0.00) | (0.00) | (0.01) | (0.00) | (0.05) | (2.98) | (0.00) | (0.01)
0.30 | 0.00 | 0.00 | 0.06 | 0.00 | 0.11 | 4.48 | 1.00 | 0.94
(0.07) | (0.00) | (0.00) | (0.02) | (0.00) | (0.07) | (6.12) | (0.00) | (0.03)
0.60 | 0.00 | 0.00 | 0.08 | 0.00 | 0.08 | 6.46 | 0.98 | 0.94
(0.05) | (0.00) | (0.00) | (0.02) | (0.00) | (0.12) | (14.51) | (0.05) | (0.03)
0.80 | 0.01 | 0.00 | 0.21 | 0.00 | 0.09 | 10.00 | 0.96 | 0.92
(0.04) | (0.00) | (0.00) | (0.04) | (0.00) | (0.08) | (21.00) | (0.06) | (0.03)
5000 | 0.15 | 0.00 | 0.00 | 0.01 | 0.00 | 0.06 | 1.89 | 0.99 | 0.93
(0.05) | (0.00) | (0.00) | (0.00) | (0.00) | (0.08) | (3.07) | (0.04) | (0.02)
0.30 | 0.00 | 0.00 | 0.03 | 0.00 | 0.09 | 5.07 | 0.94 | 0.93
(0.05) | (0.00) | (0.00) | (0.01) | (0.00) | (0.08) | (9.23) | (0.11) | (0.02)
0.60 | 0.00 | 0.00 | 0.04 | 0.00 | 0.04 | 2.05 | 0.99 | 0.94
(0.07) | (0.00) | (0.00) | (0.01) | (0.00) | (0.08) | (4.86) | (0.04) | (0.02)
0.80 | 0.00 | 0.00 | 0.10 | 0.00 | 0.10 | 12.86 | 0.98 | 0.94
(0.05) | (0.00) | (0.00) | (0.03) | (0.00) | (0.08) | (27.54) | (0.05) | (0.02)
10000 | 0.15 | 0.00 | 0.00 | 0.00 | 0.00 | 0.13 | 1.30 | 0.99 | 0.95
(0.05) | (0.00) | (0.00) | (0.00) | (0.00) | (0.21) | (2.72) | (0.04) | (0.01)
0.30 | 0.00 | 0.00 | 0.01 | 0.00 | 0.10 | 2.70 | 0.99 | 0.94
(0.08) | (0.00) | (0.00) | (0.00) | (0.00) | (0.08) | (3.88) | (0.04) | (0.03)
0.60 | 0.00 | 0.00 | 0.02 | 0.00 | 0.07 | 4.16 | 0.95 | 0.95
(0.06) | (0.00) | (0.00) | (0.00) | (0.00) | (0.09) | (9.62) | (0.06) | (0.02)
0.80 | 0.00 | 0.00 | 0.06 | 0.00 | 0.05 | 1.34 | 0.99 | 0.93
(0.04) | (0.00) | (0.00) | (0.02) | (0.00) | (0.04) | (1.61) | (0.04) | (0.03)

Table 9: Results for Pattern 1, with number of locations, $L=100$.

$N$ | Prop. of 0s | MSE | CP($\theta$)
---|---|---|---
$\beta$ | $\gamma$ | $\mathbf{w}$ | $\xi$ | $\mbox{\boldmath$\theta$}_{m}$ | $\mbox{\boldmath$\theta$}_{pr}$
2000 | 0.15 | 0.00 | 0.00 | 0.06 | 0.00 | 0.98 | 0.93
(0.01) | (0.00) | (0.00) | (0.03) | (0.00) | (0.05) | (0.03)
0.30 | 0.00 | 0.00 | 0.34 | 0.00 | 0.96 | 0.95
(0.01) | (0.00) | (0.00) | (0.10) | (0.00) | (0.06) | (0.02)
0.60 | 0.00 | 0.00 | 0.92 | 0.00 | 0.96 | 0.94
(0.01) | (0.00) | (0.00) | (0.31) | (0.00) | (0.08) | (0.03)
0.80 | 0.00 | 0.00 | 2.74 | 0.00 | 1.00 | 0.91
(0.01) | (0.00) | (0.00) | (0.26) | (0.00) | (0.00) | (0.02)
5000 | 0.15 | 0.00 | 0.00 | 0.03 | 0.00 | 0.95 | 0.95
(0.01) | (0.00) | (0.00) | (0.02) | (0.00) | (0.11) | (0.02)
0.30 | 0.00 | 0.00 | 0.13 | 0.00 | 0.98 | 0.93
(0.01) | (0.00) | (0.00) | (0.04) | (0.00) | (0.05) | (0.03)
0.60 | 0.00 | 0.00 | 0.46 | 0.00 | 1.00 | 0.95
(0.01) | (0.00) | (0.00) | (0.15) | (0.00) | (0.00) | (0.02)
0.80 | 0.00 | 0.00 | 1.68 | 0.00 | 0.91 | 0.91
(0.01) | (0.00) | (0.00) | (0.28) | (0.00) | (0.10) | (0.03)
10000 | 0.155 | 0.00 | 0.00 | 0.01 | 0.00 | 0.92 | 0.93
(0.01) | (0.00) | (0.00) | (0.01) | (0.00) | (0.06) | (0.03)
0.30 | 0.00 | 0.00 | 0.09 | 0.00 | 0.99 | 0.93
(0.01) | (0.00) | (0.00) | (0.03) | (0.00) | (0.04) | (0.03)
0.60 | 0.00 | 0.00 | 0.30 | 0.00 | 0.96 | 0.93
(0.01) | (0.00) | (0.00) | (0.07) | (0.00) | (0.06) | (0.02)
0.80 | 0.00 | 0.00 | 1.37 | 0.00 | 0.98 | 0.93
(0.01) | (0.00) | (0.00) | (0.19) | (0.00) | (0.05) | (0.03)

Table 10: Results for Pattern 2, with number of locations, $L=100$.

###  8.2 Convergence

To assess convergence we resort to monitoring trace plots, posterior density
and auto-correlation (ACF) plots for model parameters and coefficients. We
show them particularly for the index parameter $\xi$ in Figure 5 below. They
were produced using parameter settings outlined in the examples within the
R-package. Additional plots for $\beta$, $\gamma$, $\sigma^{2}$ and $\phi_{s}$
are included in examples within the R-package described in the next
subsection.

Figure 5: Convergence diagnostics for posterior samples of the index
parameter, $\xi$. Each row corresponds to a model—(first row) M1 (second row)
M2 (third row) M3 (fourth row) M4. In each row we show (left) trace, (center)
posterior probability density and (right) ACF plots.

###  8.3 R-package: sptwdglm

Algorithms featuring posteriors in Section 7 constitute our implementation
that is available in the publicly accessible open-source repository:
https://github.com/arh926/sptwdglm. It is composed of four Markov chain Monte
Carlo samplers as described below. Each function is accompanied with examples
in the package.

  1. 1.

dglm.autograd.R: Fits a DGLM to the data. It requires the response as a
vector, $y$, covariates as matrices $\mathbf{X}$ and $\mathbf{Z}$, upper and
lower bounds for the index parameter $\xi$ as input at the minimum to produce
an output. Other optional parameters are also provided for custom fine-tuning
which have optimized defaults for ease of use. This is model M1 in the
manuscript.

  2. 2.

ssdglm.autograd.R: Fits a DGLM featuring variable selection via the spike and
slab prior to the data. It requires the response as a vector, $y$, covariates
as matrices $\mathbf{X}$ and $\mathbf{Z}$, upper and lower bounds for the
index parameter $\xi$ as input at the minimum to produce an output. Other
optional parameters are also provided for custom fine-tuning which have
optimized defaults for ease of use. This is model M2 in the manuscript.

  3. 3.

spdglm.autograd.R: Fits a spatial DGLM to the data. It requires the
coordinates as a matrix, the response as a vector, $y$, covariates as matrices
$\mathbf{X}$ and $\mathbf{Z}$, upper and lower bounds for the index parameter
$\xi$ as input at the minimum to produce an output. Other optional parameters
are also provided for custom fine-tuning which have optimized defaults for
ease of use. This is model M3 in the manuscript.

  4. 4.

spssdglm.autograd.R: Fits a spatial DGLM featuring variable selection via the
spike and slab prior to the data. It requires the coordinates as a matrix, the
response as a vector, $y$, covariates as matrices $\mathbf{X}$ and
$\mathbf{Z}$, upper and lower bounds for the index parameter $\xi$ as input at
the minimum to produce an output. Other optional parameters are also provided
for custom fine-tuning which have optimized defaults for ease of use. This is
model M4 in the manuscript.

  5. 5.

FDR.R: Computes local false discovery rates (see page 10 of manuscript) for
posterior samples of $\beta$ and $\gamma$ arising from M2 or M4. Optional
parameters to be specified are threshold defaulting at 0.05 and percentage at
which FDR is computed defaulting to 5%.

##  9 Automobile Insurance Premiums, Connecticut, 2008: Additional results

Parameters | Levels | Median | Mean | SD | Lower HPD | Upper HPD
---|---|---|---|---|---|---
(Intercept) | – | 5.039 | 4.994 | 0.372 | 4.286 | 5.676
age.car | 0 | -0.373 | -0.344 | 0.224 | -0.787 | 0.111
1 | -0.306 | -0.271 | 0.218 | -0.641 | 0.224
2 | -0.196 | -0.164 | 0.219 | -0.519 | 0.333
3 | -0.371 | -0.331 | 0.262 | -0.740 | 0.262
4 | -0.187 | -0.156 | 0.225 | -0.528 | 0.369
5 | -0.111 | -0.078 | 0.227 | -0.454 | 0.437
6 | -0.350 | -0.320 | 0.230 | -0.696 | 0.213
7 | -0.401 | -0.367 | 0.252 | -0.855 | 0.138
risk | S | -0.176 | -0.173 | 0.060 | -0.292 | -0.054
agec | 2 | -0.014 | -0.013 | 0.097 | -0.208 | 0.163
3 | -0.292 | -0.292 | 0.101 | -0.504 | -0.105
4 | -0.150 | -0.152 | 0.114 | -0.347 | 0.096
5 | -0.430 | -0.424 | 0.230 | -0.881 | -0.022
gender | F | 0.451 | 0.436 | 0.220 | 0.015 | 0.838
M | 0.279 | 0.268 | 0.219 | -0.141 | 0.668
marital | M | -4.096 | -3.941 | 1.031 | -5.532 | -1.829
S | 4.954 | 4.893 | 1.751 | 2.285 | 7.588
deductible | B | 2.613 | 2.595 | 0.832 | 0.987 | 4.065
C | 1.025 | 0.999 | 0.439 | 0.049 | 1.779
D | 1.519 | 1.487 | 0.432 | 0.605 | 2.322
E | 1.626 | 1.579 | 0.429 | 0.651 | 2.375
F | 1.691 | 1.678 | 0.424 | 0.695 | 2.412
G | 1.733 | 1.717 | 0.433 | 0.879 | 2.626
H | 1.804 | 1.754 | 0.630 | 0.552 | 2.885
genderMarital | A | -4.998 | -4.940 | 1.725 | -7.626 | -2.515
B | 3.756 | 3.622 | 1.057 | 1.413 | 5.219
C | -4.521 | -4.532 | 1.823 | -7.360 | -1.766
D | 3.850 | 3.685 | 0.986 | 1.814 | 5.412

Table 11: Model coefficients for fixed effects in the mean model for M1. We
show the median, mean, standard deviation and HPDs.

Parameters | Levels | Median | Mean | SD | Lower HPD | Upper HPD
---|---|---|---|---|---|---
(Intercept) | – | 5.871 | 5.867 | 0.052 | 5.784 | 5.948
age.car | 0 | -0.814 | -0.804 | 0.088 | -0.915 | -0.663
1 | -1.494 | -1.471 | 0.136 | -1.657 | -1.217
2 | -1.122 | -1.112 | 0.118 | -1.299 | -0.918
3 | 0.720 | 0.737 | 0.162 | 0.485 | 1.009
4 | -0.827 | -0.826 | 0.072 | -0.937 | -0.709
5 | -0.045 | -0.033 | 0.122 | -0.203 | 0.170
6 | -0.502 | -0.489 | 0.057 | -0.571 | -0.392
7 | 0.672 | 0.673 | 0.128 | 0.485 | 0.887
risk | S | 1.688 | 1.693 | 0.048 | 1.620 | 1.772
agec | 2 | 1.333 | 1.343 | 0.025 | 1.317 | 1.400
3 | 1.384 | 1.388 | 0.010 | 1.374 | 1.409
4 | 1.394 | 1.396 | 0.013 | 1.370 | 1.417
5 | 1.157 | 1.160 | 0.124 | 0.981 | 1.351
gender | F | 0.168 | 0.151 | 0.053 | 0.055 | 0.214
M | -0.074 | -0.075 | 0.027 | -0.118 | -0.024
marital | M | 0.022 | 0.044 | 0.055 | -0.034 | 0.131
S | 2.308 | 2.277 | 0.243 | 1.829 | 2.665
deductible | B | -0.301 | -0.280 | 0.052 | -0.349 | -0.184
C | -0.122 | -0.131 | 0.093 | -0.306 | -0.011
D | -0.519 | -0.509 | 0.154 | -0.746 | -0.241
E | -0.297 | -0.287 | 0.183 | -0.563 | 0.017
F | 0.446 | 0.444 | 0.051 | 0.358 | 0.515
G | 1.664 | 1.662 | 0.039 | 1.596 | 1.720
H | 0.513 | 0.516 | 0.042 | 0.437 | 0.583
genderMarital | A | -0.912 | -0.904 | 0.166 | -1.175 | -0.624
B | 0.270 | 0.270 | 0.005 | 0.257 | 0.279
C | -0.660 | -0.654 | 0.126 | -0.856 | -0.456
D | -0.222 | -0.236 | 0.042 | -0.309 | -0.176
$\xi$ | – | 1.674 | 1.674 | 0.007 | 1.664 | 1.687

Table 12: Model coefficients for fixed effects in the dispersion model for M1.
We show the median, mean, standard deviation and HPDs.

Parameters | Levels | Median | Mean | SD | Lower HPD | Upper HPD
---|---|---|---|---|---|---
(Intercept) | – | 5.325 | 5.343 | 0.102 | 5.173 | 5.542
age.car | 0 | 0.023 | 0.021 | 0.080 | -0.109 | 0.163
1 | 0.084 | 0.081 | 0.080 | -0.066 | 0.215
2 | 0.194 | 0.192 | 0.080 | 0.051 | 0.330
3 | -0.036 | -0.040 | 0.080 | -0.187 | 0.094
4 | 0.021 | 0.020 | 0.080 | -0.126 | 0.155
5 | 0.072 | 0.068 | 0.080 | -0.074 | 0.209
6 | -0.022 | -0.025 | 0.081 | -0.168 | 0.116
7 | -0.072 | -0.078 | 0.083 | -0.224 | 0.057
risk | S | -0.195 | -0.195 | 0.024 | -0.238 | -0.143
agec | 2 | -0.135 | -0.135 | 0.032 | -0.200 | -0.075
3 | -0.479 | -0.480 | 0.031 | -0.540 | -0.421
4 | -0.626 | -0.626 | 0.036 | -0.696 | -0.559
5 | -0.694 | -0.694 | 0.059 | -0.814 | -0.581
gender | F | 0.602 | 0.605 | 0.090 | 0.438 | 0.789
M | 0.721 | 0.726 | 0.094 | 0.552 | 0.900
marital | M | 1.527 | 1.623 | 0.286 | 1.197 | 2.152
S | -0.354 | -0.339 | 0.258 | -0.763 | 0.222
deductible | B | 1.371 | 1.366 | 0.222 | 0.921 | 1.815
C | 0.147 | 0.088 | 0.253 | -0.353 | 0.495
D | 0.535 | 0.465 | 0.247 | 0.010 | 0.829
E | 0.718 | 0.637 | 0.243 | 0.207 | 1.010
F | 0.680 | 0.598 | 0.243 | 0.178 | 0.979
G | 0.680 | 0.599 | 0.240 | 0.188 | 0.970
H | 0.766 | 0.791 | 0.331 | 0.143 | 1.436
genderMarital | A | 0.226 | 0.215 | 0.223 | -0.240 | 0.623
B | -2.005 | -2.065 | 0.257 | -2.551 | -1.664
C | 0.281 | 0.275 | 0.273 | -0.278 | 0.720
D | -2.023 | -2.107 | 0.299 | -2.604 | -1.680

Table 13: Model coefficients for fixed effects in the mean model for M3.

Parameters | Levels | Median | Mean | SD | Lower HPD | Upper HPD
---|---|---|---|---|---|---
(Intercept) | – | 6.655 | 6.645 | 0.093 | 6.473 | 6.842
age.car | 0 | -0.837 | -0.838 | 0.049 | -0.926 | -0.734
1 | -1.053 | -1.051 | 0.044 | -1.139 | -0.967
2 | -0.916 | -0.915 | 0.042 | -0.992 | -0.828
3 | -0.889 | -0.889 | 0.041 | -0.966 | -0.805
4 | -0.875 | -0.873 | 0.041 | -0.955 | -0.793
5 | -0.824 | -0.824 | 0.040 | -0.908 | -0.750
6 | -0.798 | -0.796 | 0.040 | -0.869 | -0.712
7 | -0.764 | -0.763 | 0.040 | -0.851 | -0.693
risk | S | 0.105 | 0.104 | 0.015 | 0.075 | 0.134
agec | 2 | -0.027 | -0.027 | 0.022 | -0.067 | 0.018
3 | -0.042 | -0.042 | 0.014 | -0.072 | -0.015
4 | -0.078 | -0.078 | 0.014 | -0.105 | -0.049
5 | -0.280 | -0.279 | 0.029 | -0.338 | -0.225
gender | F | 0.004 | 0.003 | 0.038 | -0.071 | 0.077
M | 0.005 | 0.005 | 0.043 | -0.076 | 0.089
marital | M | -3.145 | -3.154 | 0.180 | -3.460 | -2.839
S | -3.280 | -3.281 | 0.205 | -3.642 | -2.900
deductible | B | -1.226 | -1.224 | 0.133 | -1.455 | -0.943
C | -0.751 | -0.758 | 0.082 | -0.940 | -0.620
D | -0.808 | -0.815 | 0.080 | -0.973 | -0.664
E | -0.753 | -0.759 | 0.078 | -0.926 | -0.626
F | -0.598 | -0.604 | 0.078 | -0.774 | -0.473
G | -0.236 | -0.240 | 0.076 | -0.405 | -0.107
H | 0.224 | 0.224 | 0.138 | -0.049 | 0.473
genderMarital | A | 3.172 | 3.169 | 0.202 | 2.800 | 3.536
B | 3.233 | 3.247 | 0.175 | 2.948 | 3.540
C | 3.307 | 3.305 | 0.203 | 2.909 | 3.639
D | 3.265 | 3.280 | 0.176 | 2.980 | 3.575
$\xi$ | – | 1.671 | 1.671 | 0.004 | 1.664 | 1.679

Table 14: Model coefficients for fixed effects in the dispersion model for M3.

## References

  * Abramowitz et al., (1988) Abramowitz, M., Stegun, I. A., and Romer, R. H. (1988).  Handbook of mathematical functions with formulas, graphs, and mathematical tables.
  * Agarwal et al., (2002) Agarwal, D. K., Gelfand, A. E., and Citron-Pousty, S. (2002).  Zero-inflated models with application to spatial count data.  Environmental and Ecological statistics, 9(4):341–355.
  * Akaike, (1998) Akaike, H. (1998).  Information theory and an extension of the maximum likelihood principle.  In Selected papers of hirotugu akaike, pages 199–213. Springer.
  * Banerjee and Carlin, (2004) Banerjee, S. and Carlin, B. P. (2004).  Parametric spatial cure rate models for interval-censored time-to-relapse data.  Biometrics, 60(1):268–275.
  * Banerjee et al., (2014) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2014).  Hierarchical modeling and analysis for spatial data.
  * (6) Berger, J. O., De Oliveira, V., and Sansó, B. (2001a).  Objective Bayesian analysis of spatially correlated data.  Journal of the American Statistical Association, 96(456):1361–1374.
  * (7) Berger, J. O., Pericchi, L. R., Ghosh, J., Samanta, T., De Santis, F., Berger, J., and Pericchi, L. (2001b).  Objective Bayesian methods for model selection: Introduction and comparison.  Lecture Notes-Monograph Series, pages 135–207.
  * Berliner, (2000) Berliner, M. (2000).  Hierarchical Bayesian modeling in the environmental sciences.  AStA Advances in Statistical Analysis, 2(84):141–153.
  * Best et al., (2000) Best, N. G., Ickstadt, K., and Wolpert, R. L. (2000).  Spatial Poisson regression for health and exposure data measured at disparate resolutions.  Journal of the American statistical association, 95(452):1076–1088.
  * Bradley et al., (2018) Bradley, J. R., Holan, S. H., and Wikle, C. K. (2018).  Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data (with discussion).  Bayesian Analysis, 13(1):253–310.
  * Bradley et al., (2020) Bradley, J. R., Holan, S. H., and Wikle, C. K. (2020).  Bayesian hierarchical models with conjugate full-conditional distributions for dependent data from the natural exponential family.  Journal of the American Statistical Association, 115(532):2037–2052.
  * Carlin and Louis, (2008) Carlin, B. P. and Louis, T. A. (2008).  Bayesian methods for data analysis.  CRC press.
  * Carvalho et al., (2010) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2010).  The horseshoe estimator for sparse signals.  Biometrika, 97(2):465–480.
  * Clark and Gelfand, (2006) Clark, J. and Gelfand, A. (2006).  A future for models and data in ecology.  Trends in Ecology and Evolution, 21:375–380.
  * Cressie, (2015) Cressie, N. (2015).  Statistics for spatial data.  John Wiley & Sons.
  * Dey et al., (2000) Dey, D. K., Ghosh, S. K., and Mallick, B. K. (2000).  Generalized linear models: A Bayesian perspective.  CRC Press.
  * Diggle et al., (1998) Diggle, P. J., Tawn, J. A., and Moyeed, R. A. (1998).  Model-based geostatistics.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 47(3):299–350.
  * Dunn and Smyth, (2005) Dunn, P. K. and Smyth, G. K. (2005).  Series evaluation of Tweedie exponential dispersion model densities.  Statistics and Computing, 15(4):267–280.
  * Dunn and Smyth, (2008) Dunn, P. K. and Smyth, G. K. (2008).  Evaluation of Tweedie exponential dispersion model densities by Fourier inversion.  Statistics and Computing, 18(1):73–86.
  * Eidsvik et al., (2012) Eidsvik, J., Finley, A. O., Banerjee, S., and Rue, H. (2012).  Approximate Bayesian inference for large spatial datasets using predictive process models.  Computational Statistics & Data Analysis, 56(6):1362–1380.
  * Fan and Liao, (2014) Fan, J. and Liao, Y. (2014).  Endogeneity in high dimensions.  Annals of statistics, 42(3):872.
  * Finley et al., (2009) Finley, A. O., Banerjee, S., and McRoberts, R. E. (2009).  Hierarchical spatial models for predicting tree species assemblages across large domains.  The annals of applied statistics, 3(3):1052.
  * Gelfand et al., (1996) Gelfand, A., Sahu, S., and Carlin, B. (1996).  Efficient Parametrizations for Generalized Linear Mixed Models.  Bayesian Statistics 5: Proceedings of the Fifth Valencia International Meeting, pages 165–180.
  * Gelfand, (2000) Gelfand, A. E. (2000).  Modeling and inference for point-referenced binary spatial data.  Generalized linear models: a Bayesian perspective, page 373.
  * Gelfand and Banerjee, (2017) Gelfand, A. E. and Banerjee, S. (2017).  Bayesian modeling and analysis of geostatistical data.  Annual review of statistics and its application, 4:245–266.
  * Gelfand et al., (1995) Gelfand, A. E., Sahu, S. K., and Carlin, B. P. (1995).  Efficient parametrisations for normal linear mixed models.  Biometrika, 82(3):479–488.
  * Gelfand et al., (2005) Gelfand, A. E., Schmidt, A. M., Wu, S., Silander Jr, J. A., Latimer, A., and Rebelo, A. G. (2005).  Modelling species diversity through species level hierarchical modelling.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 54(1):1–20.
  * George and McCulloch, (1993) George, E. I. and McCulloch, R. E. (1993).  Variable selection via Gibbs sampling.  Journal of the American Statistical Association, 88(423):881–889.
  * Girolami and Calderhead, (2011) Girolami, M. and Calderhead, B. (2011).  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214.
  * Halder et al., (2021) Halder, A., Mohammed, S., Chen, K., and Dey, D. K. (2021).  Spatial Tweedie exponential dispersion models: an application to insurance rate-making.  Scandinavian Actuarial Journal, 2021(10):1017–1036.
  * Halder et al., (2022) Halder, A., Mohammed, S., Chen, K., and Dey, D. K. (2022).  Spatial risk estimation in Tweedie double generalized linear models.  Proceedings of International E-Conference on Mathematical and Statistical Sciences: A Selcuk Meeting, 2022:62.
  * Heaton et al., (2019) Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R. B., Hammerling, D., Katzfuss, M., et al. (2019).  A Case Study Competition Among Methods For Analyzing Large Spatial Data.  Journal of Agricultural, Biological and Environmental Statistics, 24(3):398–425.
  * Hoeting et al., (1999) Hoeting, J. A., Madigan, D., Raftery, A. E., and Volinsky, C. T. (1999).  Bayesian model averaging: a tutorial (with comments by M. Clyde, David Draper and EI george, and a rejoinder by the authors.  Statistical science, 14(4):382–417.
  * Ishwaran and Rao, (2005) Ishwaran, H. and Rao, J. S. (2005).  Spike and slab variable selection: frequentist and Bayesian strategies.  The Annals of Statistics, 33(2):730–773.
  * Jørgensen, (1986) Jørgensen, B. (1986).  Some properties of exponential dispersion models.  Scandinavian Journal of Statistics, pages 187–197.
  * Jørgensen, (1987) Jørgensen, B. (1987).  Exponential dispersion models.  Journal of the Royal Statistical Society: Series B (Methodological), 49(2):127–145.
  * Jørgensen, (1992) Jørgensen, B. (1992).  Exponential dispersion models and extensions: A review.  International Statistical Review/Revue Internationale de Statistique, pages 5–20.
  * Jorgensen, (1997) Jorgensen, B. (1997).  The theory of dispersion models.  CRC Press.
  * Kokonendji et al., (2021) Kokonendji, C. C., Bonat, W. H., and Abid, R. (2021).  Tweedie regression models and its geometric sums for (semi‐) continuous data.  Wiley Interdisciplinary Reviews: Computational Statistics, 13(1):1496.
  * Lawson, (2018) Lawson, A. B. (2018).  Bayesian disease mapping: hierarchical modeling in spatial epidemiology.  Chapman and Hall/CRC.
  * Lee and Nelder, (2006) Lee, Y. and Nelder, J. A. (2006).  Double hierarchical generalized linear models (with discussion).  Journal of the Royal Statistical Society: Series C (Applied Statistics), 55(2):139–185.
  * Li and Lin, (2010) Li, Q. and Lin, N. (2010).  The Bayesian elastic net.  Bayesian analysis, 5(1):151–170.
  * Liang et al., (2008) Liang, F., Paulo, R., Molina, G., Clyde, M. A., and Berger, J. O. (2008).  Mixtures of g-priors for Bayesian variable selection.  Journal of the American Statistical Association, 103(481):410–423.
  * Mallick et al., (2022) Mallick, H., Chatterjee, S., Chowdhury, S., Chatterjee, S., Rahnavard, A., and Hicks, S. C. (2022).  Differential expression of single-cell RNA-seq data using Tweedie models.  Statistics in medicine, 41(18):3492–3510.
  * Martino et al., (2011) Martino, S., Akerkar, R., and Rue, H. (2011).  Approximate Bayesian inference for survival models.  Scandinavian Journal of Statistics, 38(3):514–528.
  * Matérn, (2013) Matérn, B. (2013).  Spatial variation, volume 36.  Springer Science & Business Media.
  * Mitchell and Beauchamp, (1988) Mitchell, T. J. and Beauchamp, J. J. (1988).  Bayesian variable selection in linear regression.  Journal of the american statistical association, 83(404):1023–1032.
  * Mohammed et al., (2021) Mohammed, S., Bharath, K., Kurtek, S., Rao, A., and Baladandayuthapani, V. (2021).  RADIOHEAD: Radiogenomic analysis incorporating tumor heterogeneity in imaging through densities.  The Annals of Applied Statistics, 15(4):1808–1830.
  * Morris et al., (2008) Morris, J. S., Brown, P. J., Herrick, R. C., Baggerly, K. A., and Coombes, K. R. (2008).  Bayesian analysis of mass spectrometry proteomic data using wavelet-based functional mixed models.  Biometrics, 64(2):479–489.
  * Nelder and Pregibon, (1987) Nelder, J. A. and Pregibon, D. (1987).  An extended quasi-likelihood function.  Biometrika, 74(2):221–232.
  * Park and Casella, (2008) Park, T. and Casella, G. (2008).  The Bayesian lasso.  Journal of the American Statistical Association, 103(482):681–686.
  * Raftery et al., (1997) Raftery, A. E., Madigan, D., and Hoeting, J. A. (1997).  Bayesian model averaging for linear regression models.  Journal of the American Statistical Association, 92(437):179–191.
  * Roberts and Stramer, (2002) Roberts, G. O. and Stramer, O. (2002).  Langevin diffusions and Metropolis-Hastings algorithms.  Methodology and computing in applied probability, 4(4):337–357.
  * Shono, (2008) Shono, H. (2008).  Application of the Tweedie distribution to zero-catch data in CPUE analysis.  Fisheries Research, 93(1-2):154–162.
  * Smyth, (1989) Smyth, G. K. (1989).  Generalized linear models with varying dispersion.  Journal of the Royal Statistical Society: Series B (Methodological), 51(1):47–60.
  * Smyth and Jørgensen, (2002) Smyth, G. K. and Jørgensen, B. (2002).  Fitting Tweedie’s compound Poisson model to insurance claims data: dispersion modelling.  ASTIN Bulletin: The Journal of the IAA, 32(1):143–157.
  * Smyth and Verbyla, (1999) Smyth, G. K. and Verbyla, A. P. (1999).  Adjusted likelihood methods for modelling dispersion in generalized linear models.  Environmetrics: The official journal of the International Environmetrics Society, 10(6):695–709.
  * Swallow et al., (2016) Swallow, B., Buckland, S. T., King, R., and Toms, M. P. (2016).  Bayesian hierarchical modelling of continuous non-negative longitudinal data with a spike at zero: An application to a study of birds visiting gardens in winter.  Biometrical Journal, 58(2):357–371.
  * Tweedie et al., (1984) Tweedie, M. C. et al. (1984).  An index which distinguishes between some important exponential families.  In Statistics: Applications and new directions: Proc. Indian statistical institute golden Jubilee International conference, volume 579, page 6o4.
  * Verbyla, (1993) Verbyla, A. P. (1993).  Modelling variance heterogeneity: residual maximum likelihood and diagnostics.  Journal of the Royal Statistical Society: Series B (Methodological), 55(2):493–508.
  * Williams and Rasmussen, (2006) Williams, C. K. and Rasmussen, C. E. (2006).  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA.
  * Wolpert and Ickstadt, (1998) Wolpert, R. L. and Ickstadt, K. (1998).  Poisson/gamma random field models for spatial statistics.  Biometrika, 85(2):251–267.
  * Yang et al., (2018) Yang, Y., Qian, W., and Zou, H. (2018).  Insurance premium prediction via gradient tree-boosted Tweedie compound poisson models.  Journal of Business & Economic Statistics, 36(3):456–470.
  * Ye et al., (2021) Ye, T., Lachos, V. H., Wang, X., and Dey, D. K. (2021).  Comparisons of zero-augmented continuous regression models from a bayesian perspective.  Statistics in Medicine, 40(5):1073–1100.
  * Zeger and Karim, (1991) Zeger, S. L. and Karim, M. R. (1991).  Generalized linear models with random effects; a Gibbs sampling approach.  Journal of the American statistical association, 86(413):79–86.
  * Zhang, (2002) Zhang, H. (2002).  On estimation and prediction for spatial generalized linear mixed models.  Biometrics, 58(1):129–136.
  * Zhang, (2013) Zhang, Y. (2013).  Likelihood-based and Bayesian methods for Tweedie compound Poisson linear mixed models.  Statistics and Computing, 23(6):743–757.
  * Zhou and Hanson, (2015) Zhou, H. and Hanson, T. (2015).  Bayesian spatial survival models.  Nonparametric Bayesian Inference in Biostatistics, pages 215–246.