Simultaneous Transformation and Rounding (STAR) Models for Integer-Valued Data

TL;DR

The paper introduces the STAR framework, a flexible Bayesian modeling approach that transforms and rounds continuous data to effectively model various types of integer-valued data, including zero-inflation and overdispersion.

Contribution

It develops a novel Bayesian framework for integer-valued data using simultaneous transformation and rounding, extending continuous data models to discrete settings.

Findings

STAR models outperform existing methods in predictive accuracy.

The Bayesian BART model effectively captures complex data distributions.

Application to healthcare and ecological data demonstrates practical utility.

Abstract

We propose a simple yet powerful framework for modeling integer-valued data, such as counts, scores, and rounded data. The data-generating process is defined by Simultaneously Transforming and Rounding (STAR) a continuous-valued process, which produces a flexible family of integer-valued distributions capable of modeling zero-inflation, bounded or censored data, and over- or underdispersion. The transformation is modeled as unknown for greater distributional flexibility, while the rounding operation ensures a coherent integer-valued data-generating process. An efficient MCMC algorithm is developed for posterior inference and provides a mechanism for adaptation of successful Bayesian models and algorithms for continuous data to the integer-valued data setting. Using the STAR framework, we design a new Bayesian Additive Regression Tree (BART) model for integer-valued data, which demonstrates impressive predictive distribution accuracy for both synthetic data and a large healthcare utilization dataset. For interpretable regression-based inference, we develop a STAR additive model, which offers greater flexibility and scalability than existing integer-valued models. The STAR additive model is applied to study the recent decline in Amazon river dolphins.

Figures

Tables

Table 1. Table 1: WAIC and seconds per 10000 effective samples for additive models (dolphins data).

	am-id	am-log	am-star-id	am-star-bc	am-star-np	am-Pois	am-NB
WAIC	2171	1964	2000	1924	1916	2821	1931
Sec / 10000 ES	32	30	52	67	95	843	295

Table 2. Table 2: WAIC for the ships , roaches , and Highway datasets.

	Model	Gauss-id	Gauss-log	star-log	star-sqrt	star-id	star-bc	star-np	Pois	NB
ships	lm	292	204	204	195	220	196	194	239	203
	am	240	188	193	173	182	175	171	-	-
	bart	242	204	200	182	191	196	187	-	-
roaches	lm	2722	1886	1791	1772	1952	1756	1759	12565	1793
	am	2711	1838	1740	1732	1928	1710	1729	-	-
	bart	2700	1833	1736	1732	1927	1708	1719	-	-
Highway	lm	495	481	481	484	496	486	479	1587	496
	am	487	465	465	473	484	478	460	-	-
	bart	466	463	465	465	469	476	464	-	-

Equations

$$y=h(y^{\ast }),$$

$$g(y^{\ast })=z^{\ast },z^{\ast }\sim {\Pi }_{\theta },$$

$$z^{\ast }(x)=\mu (x)+ϵ(x),ϵ(x)\stackrel{indep}{\sim }N(0,{\sigma }^2(x)),$$

$$\mathbb{P}\{y(x)=j\}=\mathbb{P}\left\{y^{\ast }(x)\in {\mathcal{A}}_j\right\}=\mathbb{P}\left\{z^{\ast }(x)\in g({\mathcal{A}}_j)\right\}.$$

$$\mathbb{P}\{y(x)=j\}=\Phi \left(\frac{g(a_{j+1})-\mu (x)}{\sigma (x)}\right)-\Phi \left(\frac{g(a_j)-\mu (x)}{\sigma (x)}\right).$$

$$g(t;\lambda )=\{\text{mbox}sgn(t)|t{|}^{\lambda }-1\}/\lambda ,\lambda >0$$

$$g(t)=b_I^{\prime }(t)\gamma ,$$

$${\gamma }_{\ell }={\tilde{\gamma }}_{\ell }/\sum _{k=1}^L{\tilde{\gamma }}_k,{\tilde{\gamma }}_{\ell }\stackrel{indep}{\sim }{N}_{+}({\mu }_{{\gamma }_{\ell }},{\sigma }_{\gamma }^2),\ell =1,\dots ,L,$$

$$\mu (x)=u^{\prime }\beta +\sum _j{f}_j(v_j),$$

$$\mu (x)=\sum _{k=1}^{m}f(x;T_k,M_k),$$

$$\text{mbox}{lpd}_{\mathcal{M}}(y)=\sum _{i=1}^n\log \left(\frac{1}{S}\sum _{s=1}^S{p}_{\mathcal{M}}(y_i\mid {\theta }^s)\right)$$

$$\text{mbox}{lpd}_{\text{textsc}star}(y)=\sum _{i=1}^n\log \left(\frac{1}{S}\sum _{s=1}^S\Phi \left\{\frac{g^s(a_{y_i+1})-{\mu }^s(x_i)}{{\sigma }^s(x_i)}\right\}-\Phi \left\{\frac{g^s(a_{y_i})-{\mu }^s(x_i)}{{\sigma }^s(x_i)}\right\}\right).$$

$$f(x)=10\sin (\pi x_1{x}_2)+20(x_3-0.5{)}^2+10x_4+5x_5$$

$$\text{mbox}{lpd}_{\mathcal{M}}(\tilde{y})=\frac{1}{n_{test}}\sum _{i=1}^{n_{test}}\log p_{\mathcal{M}}({\tilde{y}}_i|y)\approx \frac{1}{n_{test}}\sum _{i=1}^{n_{test}}\log \left\{\frac{1}{S}\sum _{s=1}^S{p}_{\mathcal{M}}({\tilde{y}}_i\mid {\theta }^s)\right\}$$

$$\mathbb{E}\{y(x)\}=\sum _{j=0}^{\infty }j\mathbb{P}\{y(x)=j\}\approx \sum _{j=1}^{J(x)}j\mathbb{P}\{y(x)=j\},$$

Full text

KEYWORDS: additive models, BART, count data, nonparametric regression

1 Introduction

A challenging scenario for prediction and inference occurs when the outcome variables are integer-valued, such as counts, (test) scores, or rounded data. Integer-valued data are ubiquitous in many fields, including epidemiology (Osthus et al.,, 2018; Kowal,, 2019), ecology (Dorazio et al.,, 2005), and insurance (Bening and Korolev,, 2012), among many others (Cameron and Trivedi,, 2013). Counts often serve as an indicator of demand, such as the demand for medical services (Deb and Trivedi,, 1997), emergency medical services (Matteson et al.,, 2011), and call center access (Shen and Huang,, 2008). Integer-valued data are discrete data, and exhibit a variety of complex distributional features including zero-inflation, skewness, over- or underdispersion, and in some cases may be bounded or censored. Consequently, prediction and modeling of integer-valued data—in the presence of predictors, over time intervals, and across spatial locations—remains a significant challenge.

The most widely-used models for integer-valued data build upon the Poisson distribution. However, the limitations of the Poisson distribution are well-known: the distribution is not sufficiently flexible in practice and cannot account for zero-inflation or over- and underdispersion. A common strategy is to generalize the Poisson model by introducing additional parameters, such as the quasi-Poisson (McCullagh and Nelder,, 1989), negative-binomial (Hilbe,, 2011), zero-inflated Poisson or negative-binomial (Cameron and Trivedi,, 2013; Neelon,, 2019), lognormal Poisson (Zhou et al.,, 2012), restricted generalized Poisson (Famoye,, 1993), and Conway-Maxwell Poisson models (Shmueli et al.,, 2005; Lord et al.,, 2008; Sellers and Shmueli,, 2010). A fundamental limitation of these approaches is that the additional parameters can introduce formidable challenges for estimation and computational scalability, especially in conjunction with regression, temporal, or spatial models.

In practice, however, it is exceedingly common for the discrete nature of the data to be ignored. Practitioners often log- or square-root-transform the observed integer-valued data and subsequently apply methods designed for continuous or Gaussian data. However, transformations to Gaussianity are ineffective for small counts (Warton,, 2018), while log-transformations introduce difficulties in the presence of zeros (O’Hara and Kotze,, 2010). More broadly, these approaches are not well-defined for integer-valued data: the data-generating process for a (transformed) Gaussian model cannot produce discrete data, which immediately amplifies model misspecification, limits interpretability, and undermines the reliability of inference and predictive distributions.

To address these challenges, we propose a coherent modeling framework for integer-valued data. The process is defined by simultaneously transforming and rounding (star) a continuous-valued process. First, a continuous-valued process is specified to model the dependence between (latent) variables. We focus on conditionally Gaussian regression models, but the star framework applies more broadly. Second, the latent variables are transformed for greater distributional flexibility. While the transformation may be specified in advance, such as logarithmic or square-root, we develop both parametric and nonparametric approaches to learn the transformation from the data, which improves predictive accuracy. Lastly, the transformed latent variables are filtered through a rounding operator mapping them to the (nonnegative) integers. This construction is inspired by the popular approach of transforming count data and applying Gaussian models, yet produces a mathematically consistent and well-defined integer-valued process. Importantly, we show that star processes are not merely valid integer-valued distributions, but also flexible integer-valued distributions, and can account for zero-inflation, bounded or censored data, and over- or underdispersion.

Another major benefit of star is its computational modularity: using a simple and efficient data augmentation technique, existing computational tools for Bayesian inference under continuous data models can be used for Bayesian inference under star models. As a result, star provides a cohesive framework for seamlessly adapting state-of-the-art continuous data models and algorithms to the integer-valued data setting. Using the star framework, we design—among others—a new Bayesian Additive Regression Tree (bart) model for integer-valued data. The resulting bart-star model combines the integer-valued distributional flexibility provided by star with the predictive and computational advantages inherent to bart. For synthetic data and a large healthcare utilization dataset (Section 4-6), the predictive performance of bart-star far exceeds that of competing methods which do not include both transformation and rounding in terms of out-of-sample predictive accuracy, model adequacy, and computational scalability.

We also apply star to study the recent decline in the tucuxi dolphin population, which inhabit the Amazon River. Using field survey data conducted by da Silva et al., (2018) from 1994 to 2017, we develop a star additive model for the number of observed tucuxi dolphins, which includes a smooth regression term for important predictor variables such as the year, day-of-year, and water level. The star additive model is interpretable yet flexible, and demonstrates favorable performance in model fit and computational efficiency relative to existing integer-valued models.

The remainder of the paper is organized as follows. Section 2 introduces the star framework, develops models for the unknown transformation, describes important properties, and discusses computational details for posterior inference. Section 3 provides example star models, which are applied to simulated data in Section 4 and real data in Sections 5 and 6. Section 7 concludes. Additional simulation results and empirical comparisons are in the Appendix. Methods are implemented in the R package rstar available on GitHub.

2 Simultaneously transforming and rounding

Consider a count-valued stochastic process $y\!:\mathcal{X}\to\mathcal{N}$, where $\mathcal{X}$ may correspond to predictors, times, or spatial locations and $\mathcal{N}=\{0,\ldots,\infty\}$. Although we focus on the nonnegative integers, our procedure may be trivially modified for integer-valued data and rounded data. Our goal is to construct a joint probability distribution for $y$ that simultaneously builds upon successful approaches for continuous stochastic processes (observed on $\mathbb{R}$ or $\mathbb{R}^{+}$), yet produced a flexible and well-defined distribution on $\mathcal{N}$.

To this end, we first introduce continuous-valued process $y^{*}\!:\mathcal{X}\to\mathcal{T}$, $\mathcal{T}\subseteq\mathbb{R}$ related to the observed count-valued data $y$ via

[TABLE]

where $h\!:\mathcal{T}\to\mathcal{N}$ is a rounding operator that sets $y(x)=j$ when $y^{*}(x)\in\mathcal{A}_{j}$ and $\{\mathcal{A}_{j}\}_{j=0}^{\infty}$ is a known partition of $\mathcal{T}$. For example, we may use the floor function defined by $\mathcal{A}_{j}=[a_{j},a_{j+1})=[j,j+1)$ for $j\in\mathcal{N}$; modifications are available for zero-inflated, bounded, or censored data (Section 2.2). The process $y^{*}$ operates as a continuous proxy for the observed counts $y$, which is more convenient for modeling, yet has a simple mapping to the observable data in (1). Naturally, the properties of the count-valued process $y$ will be determined by the rounding operator $h$ and the distribution of the continuous-valued process $y^{*}$.

We propose to induce a distribution on $y^{*}$ by transforming $y^{*}$ and specifying a distribution $\Pi_{\theta}$ on the transformed scale:

[TABLE]

where $g\!:\mathcal{T}\to\mathbb{R}$ is a strictly monotone function. Model (2) is inspired by the common practice of transforming count data prior to application of continuous (Gaussian) models. However, the star framework of (1)-(2) defines an integer-valued process for $y$, in which the transformation $g$ may be modeled as unknown for greater distributional flexibility.

While star is sufficiently general to incorporate any continuous family of stochastic process $\Pi_{\theta}$ for the latent $z^{*}$, an important special case of (2) is the conditionally Gaussian regression model:

[TABLE]

where $\mu(x)$ is the conditional expectation of $z^{*}(x)$ and the errors $\epsilon(x)$ are independent but possibly heteroscedastic, conditional on $x\in\mathcal{X}$. Examples of (3) include linear and additive models (Section 3.1) and bart (Section 3.2), with extensions for mixed effects models, spatio-temporal models, dynamic linear models, and factor models, among others.

Rounding of a continuous process has appeared previously in the literature (Canale and Dunson,, 2011, 2013). Our key innovation is the coupling of the transformation (2) with the regression model (3). The transformation $g$, which we model as unknown, endows the integer-valued process $y$ with greater distributional flexibility, yet leaves model (3) unchanged. This construction allows seamless integration of Bayesian models and algorithms for continuous data of the common form (3) into the integer-valued star framework, with efficient posterior inference available via a general MCMC algorithm (Section 2.3). As demonstrated extensively in the simulations and applications (Sections 4-6), models that fail to include both rounding and transformation cannot match the predictive performance of star models.

The distribution of $y$ is completely determined by the rounding operator $h$, the transformation $g$, and the distribution $\Pi_{\theta}$. Specifically, the probability mass associated to $y(x)=j$ for each integer $j\in\mathcal{N}$ is

[TABLE]

The distribution of $z^{*}$ is given by $\Pi_{\theta}$, while $g(\mathcal{A}_{j})$ is determined by the transformation $g$ and the rounding operator $h$. For model (3) and $\mathcal{A}_{j}=[a_{j},a_{j+1})$, (4) simplifies to

[TABLE]

The distribution in (5) is related to, yet distinct from, ordinal regression (McCullagh,, 1980). In ordinal regression, each term $g(a_{j})$ in (5) is replaced by an unknown latent threshold, say $\omega_{j}$, with an ordering constraint $\omega_{j}\leq\omega_{j+1}$ for all $j$. However, the latent thresholds $\omega_{j}$ are based only on the ranks of the observed data, and therefore ignore the information contained in the numeric values of the observed counts. Furthermore, since each threshold $\omega_{j}$ is unknown, ordinal regression introduces a new parameter for each unique data value, and therefore produces a heavily-parametrized model that is challenging to estimate. By comparison, star is substantially more parsimonious: if $g$ is known, no new parameters are needed, while if $g$ is unknown, only a small number of parameters are needed (see Section 2.1).

star is fundamentally different from simply rounding the predictions from a continuous data model $\Pi_{\theta}$. Post hoc rounding ignores the discrete nature of the data in model-fitting, and consequently introduces a disconnect between the fitted model and the model used for prediction. star clearly avoids this issue, and maintains the benefits of using well-known models for continuous data while producing a coherent integer-valued predictive distribution.

2.1 The transformation $g$

The transformation $g$ is a crucial component of star. When $g(t)=t$ and $z^{*}$ is a draw from a Gaussian process, star simplifies to Canale and Dunson, (2013). However, the identity transformation is suboptimal in many cases (see Sections 5 and 6). The popularity of log-linear models for count data, especially Poisson and negative-binomial models, suggests that regression effects $\mu(x)$ are often multiplicative for count data, and that the log-transformation $g(t)=\log(t)$ may be preferable for many applications. Similarly, the square-root transformation $g(t)=\sqrt{t}$ is the variance-stabilizing transformation of the Poisson distribution, and therefore is a common choice in applications of Gaussian methods to transformed count data. Empirically, the simulation studies and real data analyses in Sections 4-6 demonstrate that the transformation $g$ provides substantial improvements in modeling flexibility and accuracy relative to an untransformed approach.

When $g$ is fixed and known, the only unknowns are the parameters $\theta$ in the distribution $\Pi_{\theta}$ of the latent data $z^{*}$ in (2). Relative to the analogous transformed continuous-valued model, say $g(y)\sim\Pi_{\theta}$, the number of parameters is the same, yet star produces a coherent integer-valued process. A fixed transformation $g$ shares some characteristics with the link function of a generalized linear model (GLM) (McCullagh and Nelder,, 1989). For GLMs, the link function maps the expectation of an exponential family distribution to $\mathbb{R}$, which is modeled using a linear predictor. By comparison, star maps the continuous-valued $y^{*}$ to $\mathbb{R}$ and under (3) models $\mathbb{E}[g\{y^{*}(x)\}|x]=\mathbb{E}\{z^{*}(x)|x\}=\mu(x)$.

For general application of star, pre-specification of a transformation $g$ is restrictive. By allowing the data to inform $g$, the implied distribution for $y$ becomes more flexible, and the risk of model misspecification is lessened. For GLMs, Mallick and Gelfand, (1994) similarly relax the assumption of a known link function, adopting a nonparametric approach. For star, we require that the functions $g$ satisfy the following properties: (i) monotonicity, which preserves the ordering of the observed integers in the transformed latent space; (ii) smoothness, which provides regularization by encouraging information-sharing among nearby values; and (iii) shrinkage toward a pre-specified transformation, such as log or square-root.

A natural parametric specification for $g$ satisfying the aforementioned criteria is the (signed) Box-Cox transformation (Box and Cox,, 1964):

[TABLE]

with $g(t;\lambda=0)=\log(t)$. Box-Cox functions are a popular choice for transforming continuous data towards Gaussianity, which in the present setting is similar to (2) when $\Pi_{\theta}$ is Gaussian. Important special cases of (6) include the (shifted) identity transformation $g(t;\lambda=1)=t-1$, the (shifted and scaled) square-root transformation $g(t;\lambda=1/2)=2\sqrt{|t|}-2$, and the log-transformation. To learn the shape of the transformation, we place a prior on $\lambda$: we recommend $\lambda\sim N(1/2,1)$ truncated to $[0,3]$, which shrinks $g$ toward the (shifted and scaled) square-root transformation.

For additional flexibility, we also consider fully nonparametric specification for $g$. Consider an I-spline basis expansion (Ramsay,, 1988) for $g$:

[TABLE]

where $b_{I}$ is an $L$-dimensional vector of I-spline basis functions and $\gamma$ are the unknown basis coefficients. Since each I-spline basis function is monotone increasing, we ensure monotonicity of $g$ by restricting the elements of $\gamma$ to be positive.

We propose a prior for $\gamma$ in (7) that simultaneously enforces monotonicity, smoothness, and shrinkage toward a pre-specified transformation. However, care must be taken to ensure identifiability of the star model and retain interpretability of the parameters $\theta$ in $\Pi_{\theta}$. For model (3), arbitrary shifting and scaling of $g$ can be matched by shifting and scaling of $\mu$ and $\sigma$. The parametric transformation (6) preserves identifiability: $g(1,\lambda)=0$ for all $\lambda$ (shift constraint) and the prior on $\lambda$ is weakly informative (scale constraint). For nonparametric $g$ in (7), we resolve the identifiability issue by fixing $g(0)=0$ (shift constraint), which is satisfied automatically due to the I-spline construction, and $\lim_{t\rightarrow\infty}g(t)=1$ (scale constraint), which is enforced by constraining $\sum_{\ell=1}^{L}\gamma_{\ell}=1$. Specifically, let

[TABLE]

where $N_{+}$ is the half-normal distribution. Clearly, $\gamma_{\ell}>0$ for each $\ell$ and $\sum_{\ell=1}^{L}\gamma_{\ell}=1$, which guarantees monotonicity of $g$ and preserves identifiability of model (3). We select the prior mean $\mu_{\gamma}=(\mu_{\gamma_{1}},\ldots,\mu_{\gamma_{L}})^{\prime}$ such that $g(t)$ is a priori centered around a parametric function of interest, such as (6) with fixed $\lambda=\lambda_{0}$, and model $\sigma_{\gamma}^{2}$ with an inverse-Gamma prior to allow the data to determine the amount of shrinkage toward the parametric function of interest. Let $\bm{t}_{g}=(0,1,\ldots,a_{\max y_{i}+1}\})^{\prime}$ and let $B_{I}$ be the I-spline basis evaluated at $\bm{t}_{g}$, so $g(\bm{t}_{g})=B_{I}\gamma$. We solve $\tilde{\mu}_{\gamma}=\arg\min_{\mu_{\gamma}}||g(\bm{t}_{g};\lambda_{0})-B_{I}\mu_{\gamma}||^{2}$ subject to $\tilde{\mu}_{\gamma_{\ell}}>0$ for $\ell=1,\ldots,L$, which is a one-time cost, and normalize $\mu_{\gamma_{\ell}}=\tilde{\mu}_{\gamma_{\ell}}/\sum_{k=1}^{L}\tilde{\mu}_{\gamma_{k}}$. In the simulations and applications of Sections 4-6, we fix $\lambda_{0}=1/2$ and model $\sigma_{\gamma}^{-2}\sim\mbox{Gamma}(0.001,0.001)$.

We use quadratic I-splines with $L=2+\min\{(\mbox{\# unique }y_{i})/4,10\}$ knots, implemented using the splines2 package in R (Wang and Yan,, 2018). Boundary knots are placed at zero and $\max\{y_{i}\}$, while the $L-2$ interior knots selected using the sample quantiles of $\{y_{i}\}$ excluding zero, one, and $\max\{y_{i}\}$, with an interior knot placed at one to improve distributional flexibility near zero. Since Ramsay, (1988) use $L=3$ or $L=5$ in all monotone spline examples, a small number of knots may be adequate in many cases.

2.2 Model properties

By design, star builds upon models for continuous data, such as those in Section 3, and adapts them for integer-valued data. Yet star is not merely a mechanism for producing valid integer-valued processes: star also provides important distributional properties for modeling integer-valued data in practice. By careful selection of the rounding operator $h$ and the transformation $g$, star provides the capability to model zero-inflation, bounded or censored data, and over- or underdispersion.

In applications with count data, it is common to observe an abundance of zeros, $y(x)=0$. star can be parametrized such that zero counts occur whenever $z^{*}(x)$ is negative:

Lemma 1** (Zero-inflation).**

For any star model with $g(\mathcal{A}_{0})=(-\infty,0)$, we have (i) $y(x)=0$ if and only if $z^{*}(x)<0$ and (ii) $\mathbb{P}\{y(x)=0\}=\mathbb{P}\{z^{*}(x)\leq 0\}.$

Lemma 1 is valid for known or unknown transformations, and is easily satisfied for (6) letting $\mathcal{A}_{0}=(a_{0},a_{1})=(-\infty,1)$ for $\lambda\neq 0$ and $\mathcal{A}_{0}=(a_{0},a_{1})=(0,1)$ for $\lambda=0$. For model (3), $\Pi_{\theta}$ is conditionally Gaussian, which may place substantial prior mass on $z^{*}(x)<0$ and thus $y(x)=0$. Therefore, star has a built-in and interpretable mechanism for handling zero counts, and does not require the addition of an artificial constant to the transformation, such as $\log(y+1)$. Furthermore, dependence among zero values is implicit in the model: $\mathbb{P}\{y(x)=0,y(x^{\prime})=0\}=\mathbb{P}\{z^{*}(x)<0,z^{*}(x^{\prime})<0\}$ depends on the joint distribution of $(z^{*}(x),z^{*}(x^{\prime}))$, which is modeled by $\Pi_{\theta}$.

Another common characteristic of count-valued data is a deterministic upper bound $K$. For instance, if $y$ counts the number of days on which an event occurred in a given year, then $y(x)\in\{0,1,\ldots,K\}$ and $K=365$. star can easily incorporate this information into the distribution for $y$ as formalized in the next lemma.

Lemma 2** (Bounded observations).**

For any star model, letting $g(\mathcal{A}_{K})=[g(a_{K}),\infty)$ implies $\mathbb{P}\{y(x)\leq K\}=1$.

The boundedness constraints in Lemma 2 are compatible with any choice of (unconstrained) continuous-valued model (2) and do not require modification of the algorithms for estimation and inference in Section 2.3. Similar to the case of zero values in Lemma 1, star allows for dependence among $y$ values that attain the upper bound: $\mathbb{P}\{y(x)=K,y(x^{\prime})=K\}=\mathbb{P}\{z^{*}(x)\geq g(a_{K}),z^{*}(x^{\prime})\geq g(a_{K})\}$, which again is modeled by (2) or (3). When $K=1$, the star model (3) with $\mu(x)=x^{\prime}\beta$ and $\sigma^{2}(x)=1$ simplifies to probit regression.

Interestingly, the construction in Lemma 2 is coherent under right-censoring, which occurs when an observed count value of $K$ implies that $y(x)\geq K$. Right-censoring is common in surveys, where large values are often grouped together. The following lemma formalizes the properties of star subject to right-censoring of the observations.

Lemma 3** (Right-censoring).**

Let $\textsc{star}(h,g,\Pi_{\theta})$ denote model (1)-(2) with rounding operator $h$, transformation $g$, and latent distribution $\Pi_{\theta}$. For right-censored observations $y_{c}(x)=\min\{y(x),K\}$ with $y\sim\textsc{star}(h,g,\Pi_{\theta})$ and $y_{c}\sim\textsc{star}(h^{\prime},g,\Pi_{\theta})$ such that $h$ and $h^{\prime}$ satisfy $a_{K}=a^{\prime}_{K}$ and $\mathcal{A}^{\prime}_{K}=[a_{K},\infty)$, we have $\mathbb{P}\{y(x)\geq K\}=\mathbb{P}\{z^{*}(x)\geq g(a_{K})\}=\mathbb{P}\{y_{c}(x)=K\}$.

For right-censored data, the likelihood includes terms of the form $\mathbb{P}\{y(x)\geq K\}$ for censored observations. Lemma 3 shows that the censored likelihood terms under a star model for $y$ are equivalent to the non-censored likelihood terms $\mathbb{P}\{y_{c}(x)=K\}$ under a star model for $y_{c}=\min\{y(x),K\}$ with $\mathcal{A}^{\prime}_{K}=[a_{K},\infty)$. Remarkably, star preserves the correct right-censored likelihood for $y$ by directly modeling the observed counts $y_{c}$ and setting $\mathcal{A}_{K}^{\prime}=[a_{K},\infty)$, with no further modifications needed for the model specification or estimation procedure. By comparison, common parametric approaches for modeling count data, such as the Poisson model and its generalizations, require careful modifications of the likelihood and tailored algorithms for estimation and inference in the case of right-censoring. Naturally, a similar approach is available for left-censoring.

Lastly, we show that star processes are capable of modeling over- or underdispersion. In Figure 1, we illustrate the relationships among the expectation $\mathbb{E}[y]$, the variance $\mbox{Var}(y)$, and the probability of zeros $\mathbb{P}(y=0)$ for a star process defined by $z^{*}\sim N(\mu,\sigma^{2})$ and transformation (6) with $\lambda=1/2$. For different values of the parameters the star process exhibits different features, including overdispersion, underdispersion, and zero-inflation.

2.3 Posterior inference

We develop a general Markov Chain Monte Carlo (MCMC) algorithm for Bayesian inference under star. The hierarchical construction of star in (1)-(2) is accompanied by a computationally convenient data augmentation strategy, which we leverage to incorporate existing sampling techniques for the unknown parameters $\theta$ in (2). To emphasize the modularity of the proposed approach, we omit model-specific details for sampling $\theta$ until Section 3.

Let $\mathcal{D}=\{x_{i},y_{i}\}_{i=1}^{n}$ denote the observed pairs of points $x_{i}\in\mathcal{X}$ and integer-valued data $y_{i}=y(x_{i})$. Consider a Bayesian specification of (2) with suitable prior on $\theta$ and an algorithm $\mathbb{A}$ which draws from the posterior distribution of $\theta$ given the (continuous) data. The sampling algorithm $\mathbb{A}$ is designed for continuous data, such as (Gaussian) additive models or bart, depending on the different choices for $\Pi_{\theta}$ (see Section 3 for details). The posterior sampling algorithm for star defines a Gibbs sampler by combining a data augmentation step with algorithm $\mathbb{A}$ as follows:

Sample $\left[z^{*}(x_{i})\mid\mathcal{D},\theta\right]$ from $\Pi_{\theta}$ truncated to $z^{*}(x_{i})\in g(\mathcal{A}_{y_{i}})$ for $i=1,\dots n$; 2. 2.

Sample $\left[\theta\mid\bm{z}^{*}\right]$ using algorithm $\mathbb{A}$ conditioning on $\bm{z}^{*}=(z^{*}(x_{1}),\ldots,z^{*}(x_{n}))^{\prime}$.

In the case of model (3), the data augmentation step may be computed efficiently using a standard univariate truncated normal sampler. Specifically, for $\mathcal{A}_{j}=[a_{j},a_{j+1})$, the full conditional distribution of the latent data is $\left[z^{*}(x_{i})\mid\mathcal{D},\theta\right]\sim N(\mu(x_{i}),\sigma^{2}(x_{i}))$ truncated to $[g(a_{y_{i}}),g(a_{y_{i}+1}))$. While the process $y^{*}$ is useful for interpretability of the star model, it is not necessary for inference or sampling.

When the transformation $g$ is unknown, an additional sampling step is required. For the parametric Box-Cox case (6) this translates to sampling the parameter $\lambda$ from its full conditional posterior distribution, for which we use a slice sampler (Neal,, 2003). For the nonparametric model (7) with prior (8), we sample $\xi_{\gamma}=\log(\tilde{\gamma})$ using Metropolis-Hastings and \big{[}\sigma_{\gamma}^{-2}\mid-\big{]}\sim\mbox{Gamma}\{0.001+L/2,0.001+\sum_{\ell=1}^{L}(\tilde{\gamma}_{\ell}-\mu_{\gamma_{\ell}})^{2}\}, and set $g(t)=b_{I}^{\prime}(t)\gamma$ as defined in (7)-(8). The sampler for $\xi_{\gamma}$ uses a Gaussian random walk proposal with covariance matrix tuned using the robust adaptive Metropolis (RAM) algorithm of Vihola, (2012) during a preliminary burn-in period. Within the RAM algorithm we set a target acceptance rate of 30% with an adaptation rate of 0.75; see Vihola, (2012) for details. We adapt the proposal covariance only during the first 50% of the burn-in period, so the MCMC draws we save for inference are generated from a (non-adaptive) Metropolis-within-Gibbs sampling algorithm.

To simulate from the posterior predictive distribution $[\tilde{y}(x)|\mathcal{D}]$, we additionally sample $\left[\tilde{z}^{*}(x)\mid\theta\right]$ from $\Pi_{\theta}$ using the current draw of $\theta$ and set $\tilde{y}(x)=h\left[g^{-1}\left\{\tilde{z}^{*}(x)\right\}\right]$ for each $x$. This step is extremely simple, yet provides inference for integer-valued predictions, model-based imputation of missing data at $x\in\mathcal{X}$, and useful model diagnostics. For parametric $g$ in (6), the functions $g^{-1}(s;\lambda)$ are known, while for nonparametric $g$ we approximate $g^{-1}(s)\approx\arg\min_{t}|s-b_{I}^{\prime}(t)\gamma|$, where the minimum $t$ is computed over a grid of values.

The proposed framework for MCMC balances modularity and flexibility: it combines existing algorithms for continuous data models with a transformation to provide distributional flexibility for integer-valued data. The importance of modularity has been demonstrated recently for the negative-binomial distribution, for which Polson et al., (2013) developed a Pólya-Gamma data augmentation scheme for Gibbs sampling. This approach has allowed a variety of Gaussian models to be extended for negative-binomial data, including linear regression (Zhou et al.,, 2012), factor models (Klami,, 2015), and functional time series models (Kowal,, 2019), yet faces two important limitations: first, it is restricted to the negative-binomial distribution, and second, the resulting MCMC sampler is often inefficient (Duan et al.,, 2018). As demonstrated in Section 6 and Appendix A.2, the proposed star algorithm provides excellent MCMC efficiency, even for nonlinear versions of model (3).

3 Regression modeling with STAR

For inference and prediction of integer-valued data $y$ observed with predictors $x$, we apply the star modeling framework to develop additive (Section 3.1) and bart (Section 3.2) regression models. Each model may be combined with a known or unknown transformation, and posterior inference proceeds using the general approach from Section 2.3. The additive and bart star models are evaluated for synthetic data in Section 4 and real data in Sections 5 and 6, with additional model comparisons and diagnostics in the Appendix.

3.1 Additive models

Suppose the $p$ predictors are partitioned as $x^{\prime}=(u^{\prime},v^{\prime})$ for linear predictors $u$ and nonlinear predictors $v$. The star additive model is given by (3) with conditional mean

[TABLE]

where $f_{j}\!:\mathcal{X}_{j}\to\mathbb{R}$ is an unknown function of $v_{j}\in\mathcal{X}_{j}$. The unknown $f_{j}$ are typically modeled as smooth nonparametric functions, and may capture nonlinearities in each $v_{j}$. The star linear model is a special case of (9) with $\mu(x)=x^{\prime}\beta$. For the conditional variance of $z^{*}$ in (3), we use the conditionally conjugate prior $\sigma^{-2}\sim\mbox{Gamma}(0.001,0.001)$.

Within the star framework, we apply flexible and computationally efficient parametrizations of (9) that have been well-developed for Gaussian and exponential family models. The linear regression coefficients are assigned conditionally Gaussian priors, $\beta\sim N(0,\Sigma_{\beta})$, including the ridge prior $\Sigma_{\beta}=\sigma_{\beta}^{2}I$ and other shrinkage priors (Carvalho et al.,, 2010) as special cases. The nonlinear functions in (9) are modeled smoothly using a basis expansion $f_{j}(v_{j})=b_{j}^{\prime}(v_{j})\alpha_{j}$, where $b_{j}$ is an $L_{j}$-dimensional vector of known basis functions and $\alpha_{j}$ is a vector of unknown coefficients. We select a cubic P-spline basis with second-order difference penalty on the coefficients, which may be reparametrized such that $\alpha_{j}\sim N(0,\sigma_{\alpha_{j}}^{2}I)$ is the smoothing prior and $B_{j}^{\prime}B_{j}$ is diagonal, where $B_{j}=(b_{j}(v_{1,j}),\ldots,b_{j}(v_{n,j}))^{\prime}$ is the basis matrix (Scheipl et al.,, 2012). The nonlinear terms are constrained such that $\sum_{i=1}^{n}f_{j}(v_{i,j})=0$ for identifiability, which is enforced in the reparametrization of the basis $b_{j}$. We let $\sigma_{\alpha_{j}}^{-2}\sim\mbox{Gamma}(0.1,0.1)$, which allows the smoothness of each $f_{j}$ to be learned from the data.

For observed predictors $x_{i}^{\prime}=(u_{i}^{\prime},v_{i}^{\prime})$, let $U$ denote the matrix of linear predictors and let $\bm{f}_{j}=B_{j}\alpha_{j}$. The Gibbs sampler for the star additive model iterates the following full conditional distributions:

Sample $[z^{*}(x_{i})\mid-]\sim N(u_{i}^{\prime}\beta+\sum_{j}f_{j}(v_{i,j}),\sigma^{2})$ truncated to $z^{*}(x_{i})\in g(\mathcal{A}_{y_{i}})$; 2. 2.

Sample [\beta\mid-]\sim N\big{(}Q_{\beta}^{-1}\ell_{\beta},Q_{\beta}^{-1}\big{)} where $Q_{\beta}=\sigma^{-2}U^{\prime}U+\Sigma_{\beta}^{-1}$ and $\ell_{\beta}=\sigma^{-2}U^{\prime}(\bm{z}^{*}-\sum_{j}\bm{f}_{j})$; 3. 3.

For each $j$, sample [\alpha_{j}\mid-]\sim N\big{(}Q_{\alpha_{j}}^{-1}\ell_{\alpha_{j}},Q_{\alpha_{j}}^{-1}\big{)} where $Q_{\alpha_{j}}=\sigma^{-2}B_{j}^{\prime}B_{j}+\sigma_{\alpha_{j}}^{-2}I$ and $\ell_{\alpha_{j}}=\sigma^{-2}B_{j}^{\prime}(\bm{z}^{*}-U\beta-\sum_{k\neq j}\bm{f}_{k})$ and set $\bm{f}_{j}=B_{j}\alpha_{j}$; 4. 4.

Sample [\sigma^{-2}\mid-]\sim\mbox{Gamma}\big{(}0.001+n/2,0.001+||\bm{z}^{*}-U\beta-\sum_{j}\bm{f}_{j}||^{2}/2\big{)}; 5. 5.

For each $j$, sample [\sigma_{\alpha_{j}}^{-2}\mid-]\sim\mbox{Gamma}\big{(}0.1+L_{j}/2,0.1+\sum_{\ell=1}^{L_{j}}\alpha_{j}^{2}/2\big{)}.

The MCMC sampling algorithm for the star additive model is efficient, with empirical support provided in Section 6. The computational complexity for the nonlinear basis coefficients $\{\alpha_{j}\}$ is $\mathcal{O}(\sum_{j}L_{j})$ due to the diagonality of $B_{j}^{\prime}B_{j}$, while the linear coefficients $\beta$ also may be sampled efficiently (Rue,, 2001; Bhattacharya et al.,, 2016). Additional sampling steps for $\Sigma_{\beta}$ depend on the model specification, but are often available in closed form.

3.2 Bayesian Additive Regression Trees

While additive models are effective at capturing nonlinear marginal effects, they are often inadequate for modeling interactions among predictors. Specific pairwise or higher order interactions may be specified in advance, but including all possible interactions in an additive model requires a massive number of parameters. As a remedy, Chipman et al., (2010) proposed bart, which is a “sum-of-trees” model within a fully Bayesian framework. Tree-based regression models, such as Chipman et al., (1998), are designed to model complex interactions among predictors. Notably, bart utilizes many trees, where each individual tree is constrained via the prior to be a weak learner. As a result, bart provides the capability to capture nonlinear interactions yet features built-in mechanisms to guard against overfitting. For continuous and binary data, the predictive performance of bart is highly competitive with state-of-the-art statistical and machine learning models.

For integer-valued data, bart has been relatively underutilized. Adaptations of bart for negative-binomial data are feasible via Pólya-Gamma augmentation (Polson et al.,, 2013), similar to the probit implementation in Chipman et al., (2010) for binary data. However, this approach is limited in distributional flexibility, and the MCMC inefficiencies of Pólya-Gamma augmentation are unlikely to be ameliorated given the complexity of the (Gaussian) bart sampling algorithm. Recently, Murray, (2017) proposed a log-linear bart model for count-valued and categorical data using a parameter expansion valid for certain likelihoods in log-linear models, in particular (zero-inflated) negative-binomial and Poisson. However, extensions to more flexible count distributions may require alternative computational strategies and appropriately modified prior distributions.

Within the star framework, we parametrize the bart model (bart-star) as in Chipman et al., (2010) and specifically

[TABLE]

where $T_{k}$ is a binary tree comprised of interior splitting rules and terminal nodes and $M_{k}=\{\eta_{1,k},\ldots,\eta_{b_{k},k}\}$ is the value at each of $b_{k}$ terminal nodes for tree $T_{k}$. For a given predictor $x$, each tree $T_{k}$ in (10) assigns a value $\eta_{\ell,k}\in M_{k}$, and these values are summed across all trees $k=1,\ldots,m$. Chipman et al., (2010) propose prior distributions that constrain each $T_{k}$ to be shallow, thereby limiting the order of interactions, and constrain each $\eta_{\ell,k}$ to be small, thereby limiting the contribution of each tree. Both mechanisms guard against overfitting, and in combination produce a sum of weak learners. The joint prior distribution is specified as a prior for the tree, $p(T_{k})$, which follows Chipman et al., (1998), and a prior for the terminal values given the tree, $p(\eta_{\ell,k}|T_{k})$, which is Gaussian. A key feature of star is that, by transforming to Gaussianity, we inherit the same framework as the original bart, and therefore may directly incorporate the well-studied priors and hyperparameters from Chipman et al., (2010).

More careful consideration is required for the prior distribution of $\sigma^{2}$ in (3). Chipman et al., (2010) emphasize that an informative prior distribution is important to balance between overly aggressive and overly conservative model fits, and parametrize the prior for $\sigma^{2}$ as an inverse chi-square distribution calibrated using a data-based overestimate $\hat{\sigma}$ of $\sigma$. However, any statistics calculated from the data $y$ are likely inappropriate for star, since the transformation $g$ impacts the scale of $z^{*}$. As a remedy, we compute $\hat{\sigma}$ as the posterior median of $\sigma$ from the star linear model (9), where the transformation in the linear model is chosen to match the transformation in bart-star. Given $\hat{\sigma}$, which indeed is a data-based overestimate of $\sigma$, we adopt the default hyperparameter suggestions of Chipman et al., (2010).

For posterior inference under bart-star, we combine a data augmentation step for $z^{*}(x_{i})$ similar to Section 3.1 with a sweep from the original Chipman et al., (2010) algorithm to draw the bart parameters in (10) using $\bm{z}^{*}$ as data. The Chipman et al., (2010) bart sampler proceeds using backfitting: draws for the $k$th tree $[(T_{k},M_{k})|\bm{z}^{*},\{(T_{k^{\prime}},M_{k^{\prime}})\}_{k\neq k^{\prime}},\sigma]$ are generated using Chipman et al., (1998) and $\sigma^{2}$ is sampled from an inverse-Gamma distribution. Incorporating these sampling steps into the larger Gibbs sampler in Section 2.3 is straightforward using the dbarts package (Dorie et al.,, 2018).

4 Simulation studies

The proposed star modeling framework is evaluated using simulated data and compared to existing methods for Poisson, negative-binomial, and Gaussian data. Synthetic data $y_{i}$ for $i=1,\dots n$ and $n=100$ are simulated from a negative-binomial distribution parametrized by conditional mean $\mathbb{E}\{y_{i}(x)|\lambda_{i}^{*}\}=\lambda_{i}^{*}(x)$ and variance $\mbox{Var}\{y_{i}(x)|\lambda_{i}^{*},r^{*}\}=\lambda_{i}^{*}(x)\left\{1+\lambda_{i}^{*}(x)/r^{*}\right\}$ with dispersion parameter $r^{*}>0$. As $r^{*}$ decreases to zero, the variance increasingly dominates the mean while as $r^{*}\rightarrow\infty$, the distribution converges to a Poisson distribution with parameter $\lambda_{i}^{*}(x)$. We select $r^{*}=1$ to simulate negative-binomial data with large overdispersion and $r^{*}=1000$ to simulate approximate Poisson data. We consider linear (Section 4.1) and nonlinear (Section 4.2) parametrizations for the log-mean $\log\lambda_{i}^{*}(x)$. We emphasize that in all cases, the simulated datasets are not generated under the proposed star model: they are simulated from negative-binomial and (approximate) Poisson distributions.

Competing models are compared using the Watanabe-Akaike/widely-applicable information criteria (WAIC) (Watanabe,, 2010). WAIC estimates out-of-sample predictive accuracy using a single model fit requiring only minimal additional computations, and is asymptotically equivalent to cross-validation. The WAIC for a model $\mathcal{M}$ is defined as $\mbox{WAIC}_{\mathcal{M}}=-2\left(\mbox{lpd}_{\mathcal{M}}-d_{\mathcal{M}}\right)$, where $d_{\mathcal{M}}$ is the effective number of parameters for model ${\mathcal{M}}$ and $\mbox{lpd}_{\mathcal{M}}$ is the log-predictive pointwise density defined by

[TABLE]

for $\theta^{s}$ drawn from its posterior distribution. For star models, we simply have

[TABLE]

For the effective number of parameters, we follow the recommendation of Gelman et al., (2014) and use the sample variance of the pointwise log-likelihoods across MCMC simulations: $d_{\mathcal{M}}=\sum_{i=1}^{n}\mbox{Var}\left(\log p(y_{i}\mid\theta^{s})\right)$. The pointwise log-likelihood of star is simple and efficient to compute, and is sufficient for computing WAIC as well as other information criteria.

WAIC is used in Sections 4.1 and 4.2 for synthetic data, but exact out-of-sample metrics are provided for the real data in Section 5. Appendix A.1 also provides root mean square errors for estimating the conditional expectation $\lambda_{i}^{*}(x)$ to evaluate point estimation accuracy.

4.1 Linear mean functions

We first consider a linear log-mean, $\log\lambda_{i}^{*}(x)=\beta_{0}+\sum_{j=1}^{p}x_{i,j}\beta_{j}$, where the $p=6$ predictors are drawn independently from $x_{i,j}\sim N(0,1)$ and the coefficients are $\beta_{0}=\log(1.5)$, $\beta_{1}=\beta_{2}=\beta_{3}=\log(2.0)$, and $\beta_{4}=\beta_{5}=\beta_{6}=0$. Under this specification, the expected counts at $x_{i,j}=0$ is 1.5, while each nonzero coefficient $\beta_{j}$ for $j=1,2,3$ increases the expected counts by a factor of 2 per one unit change in each $x_{j}$.

For comparison, we consider a variety of Bayesian linear regression models. Among star models, we use the linear model (9) with prior $[\beta_{j}|\sigma_{\beta}]\sim N(0,\sigma_{\beta}^{2})$ and $[\sigma_{\beta}]\sim\mbox{Uniform}(0,10^{4})$, and the following transformations: known transformation (6) with $\lambda=0$ (lm-star-log), $\lambda=1/2$ (lm-star-sqrt), and $\lambda=1$ (lm-star-id); unknown parametric transformation (6) (lm-star-bc); and unknown nonparametric transformation (7) (lm-star-np). We also include the same linear model, but instead with the Gaussian model (3) applied directly to the counts $y$ (lm) and the log-transformed counts $\log(y+1)$ (lm-log). These models are natural competitors to star, since they incorporate the same model for $\mu(x)$ but omit the rounding step in (1) and therefore do not produce an integer-valued distribution. Lastly, we include Poisson (lm-Pois) and negative-binomial (lm-NegBin) linear regression models with a log-link, implemented using the rstanarm package (Goodrich et al.,, 2018). lm-Pois and lm-NegBin are widely used for modeling count data, and here correspond to the true data-generating process.

In Figure 2, we plot the relative WAICs across simulated datasets, defined as the ratio between the WAIC of the generic model over the WAIC for a baseline method, for which we select lm-log. Relative WAIC standardizes model performance across simulated datasets: methods with a relative WAIC less than one demonstrate improvement relative to the baseline method. The star models, particularly those with unknown transformations (lm-star-bc and lm-star-np), offer substantial improvements relative to lm-log, and are highly competitive with the true models lm-Pois and lm-NegBin. The star model improvements relative to the Gaussian models and the identity transformation lm-star-id definitively demonstrate the importance of both rounding and transformation.

4.2 Nonlinear mean functions

To evaluate bart-star, we specify a nonlinear form for the log-mean, $\log\lambda_{i}^{*}(x)=\beta_{0}+\beta_{1}\tilde{f}(x)$, where $\tilde{f}(x)$ is the centered and scaled Friedman function (Friedman,, 1991)

[TABLE]

featured in the original bart simulations (Chipman et al.,, 2010). As in Chipman et al., (2010), we select $p=10$ and simulate $x_{i,j}\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{Uniform}(0,1)$. We fix the parameters $\beta_{0}=\log(1.5)$ and $\beta_{1}=\log(5.0)$, which again corresponds to low counts with a moderate signal.

We combine the bart-star model of Section 3.2 with known transformation (6) for $\lambda=0$ (bart-star-log), $\lambda=1/2$ (bart-star-sqrt), and $\lambda=1$ (bart-star-id); unknown parametric transformation (6) (bart-star-bc); and unknown nonparametric transformation (7) (bart-star-np). For competitors, we include the Gaussian bart model (bart-id) of Chipman et al., (2010) and a Gaussian bart model on the log-transformed counts $\log(y+1)$ (bart-log). Lastly, we include the linear models lm-star-bc and lm-log from Section 4.1.

The relative WAICs are plotted in Figure 3, where again we use the log-transformed Gaussian model (bart-log) as the baseline. Both bart-star-id and bart-id are omitted as noncompetitive, and indicates the importance of an appropriate transformation. bart-star provides substantial improvements relative to bart and linear Gaussian models, with especially strong performance from the unknown transformation models (bart-star-bc and bart-star-np). Perhaps surprisingly, the star linear model lm-star-bc outperforms bart-log for negative-binomial data, despite the nonlinearity in (12). By comparison, bart-star-bc consistently outperforms lm-star-bc, which suggests that the proposed bart-star model is capable of detecting the nonlinear features in (12).

5 Predicting the demand for healthcare utilization

Individualized prediction of healthcare utilization is critical both for assessing the health risks of an individual and for monitoring the aggregate stress on the healthcare system. By providing more accurate predictive distributions of individual healthcare utilization, it is possible to obtain uncertainty quantification for various measures of individual and aggregate demand, and consequently achieve more efficient allocation of medical resources and more informed patients. To assess the predictive ability of star models for this task, we use data from the National Medical Expenditure Survey (NMES) available in the AER package in R (Kleiber and Zeileis,, 2008). Multiple measures of healthcare utilization are available, including physician office visits (visits), non-physician office visits (nvisits), physician hospital outpatient visits (ovisits), and non-physician hospital outpatient visits (novisits). Individualized predictors are also provided, including health measures, socioeconomic and demographic variables, and indicators of each patient’s type of insurance. We consider a subset of $n=4406$ elderly adults (aged 66 and older) covered by Medicare, which was previously analyzed by Deb and Trivedi, (1997) and Cameron and Trivedi, (2013).

The NMES data provides a unique opportunity for insightful out-of-sample prediction comparisons. Each measure of healthcare utilization (visits, nvisits, ovisits, and novisits) is count-valued with distinct characteristics: the probability mass functions in Figure 4 illustrate the differences in the marginal distributions, most notably the proportion of zeros and the degree of overdispersion. An adequate prediction of individual healthcare utilization may require prediction of one or more of these response variables, each of which presents unique count-valued distributional features, and which share a common set of $p=17$ individual predictor variables.

We consider out-of-sample predictive distribution accuracy for each response variable in Figure 4. In all cases, we select $n_{train}=3525$ (80%) individuals randomly for training and evaluate the predictive accuracy for the remaining $n_{test}=881$ test individuals, and repeat this exercise for 100 iterations. Posterior predictive distributions were computed for each model in Section 4.2; for conciseness, we report results for bart-log, bart-star-id, bart-star-log, and bart-star-bc. The (untransformed) bart model of Chipman et al., (2010) was noncompetitive and is omitted from the subsequent results.

To evaluate the aggregate predictive distribution accuracy, we compute the log-predictive density score for the test data $\{\tilde{y}_{i}\}_{i=1}^{n_{test}}$ for each model $\mathcal{M}$:

[TABLE]

where $\theta^{s}\sim p_{\mathcal{M}}(\theta|y)$ is a draw from the posterior under model $\mathcal{M}$. For star with model (3), we have $p_{\mathcal{M}}(\tilde{y}_{i}\mid\theta^{s})=\Phi\left[\{g^{s}(a_{\tilde{y}_{i}+1})-\mu^{s}(x_{i})]\}/\sigma^{s}(x_{i})\right]-\Phi\left[\{g^{s}(a_{\tilde{y}_{i}})-\mu^{s}(x_{i})\}/\sigma^{s}(x_{i})\right]$ similar to (11). The results for each response variable are in Figure 5; larger values indicate more accurate predictive distributions. The star model with unknown transformation (bart-star-bc) performs best in all cases, closely followed by the star model with known transformation (bart-star-log). Omitting the transformation (bart-star-id), similar to Canale and Dunson, (2013), leads to substantial deterioration in predictive accuracy. Furthermore, the untransformed bart-star-id model produced infinite scores for some test points $\tilde{y}_{i}$; these points were excluded for computing (13), but imply that bart-star-id performs even worse than indicated.

For a more targeted assessment of the predictive distributions, we compare the precision and coverage of the 90% prediction intervals for each method. Interval precision is measured by the mean prediction interval width (MPIW) computed across all $\{\tilde{y}_{i}\}_{i=1}^{n_{test}}$: smaller intervals that provide the correct coverage are preferable. The MPIWs with empirical coverages are displayed in Figure 6. Across all responses, the star model with unknown transformation (bart-star-bc) provides the most precise prediction intervals with correct coverage. For the responses with greater zero-inflation and overdispersion (nvisits, ovisits, and novisits), bart-star-bc reduces MPIWs by a median of 45%, 50%, and 59%, respectively, relative to bart-star-id, which indicates a substantial gain in predictive precision offered by the (unknown) transformation.

Lastly, we consider a specific prediction task of interest: estimating the probability that an individual will utilize the healthcare system, $\mathbb{P}(\tilde{y}_{i}>0|y)$. We evaluate each method using logarithmic scoring on the event $\{\tilde{y}_{i}>0\}$, which is a proper scoring rule for binary outcomes (Gneiting and Raftery,, 2007). The results are in Figure 7, where larger values indicate superior predictive accuracy. The bart-star models which include both transformation and rounding (bart-star-log and bart-star-bc) decisively outperform the competitors, especially the continuous-valued model (bart-log).

Based on multiple measures of predictive accuracy for each of four count-valued measures of healthcare utilization, the results from our out-of-sample comparison provide definitive confirmation of the predictive capabilities of star models, and in particular bart-star models which include both transformation and rounding.

6 Modeling the decline in Amazon river dolphins

The tucuxi dolphin (sotalia fluviatilis) is a small river dolphin that inhabits the Amazon River. While the tucuxi dolphin population was once stable, the progression of habitat degradation, dolphin fishing, and other human interference has to led to increased concerns of population decline. To assess the validity of these concerns, da Silva et al., (2018) gathered data from 1994 to 2017 using multiple observers to search for tucuxi dolphins along a particular segment of the Amazon River. In addition to the number of tucuxi dolphins observed, the data include the water level (in meters), the number of observers present, and the date for each of $n=312$ surveys. While da Silva et al., (2018) fit a linear model to the logarithm of dolphin counts, we propose to leverage the star modeling framework to investigate nonlinear effects and provide greater integer-valued distributional flexibility.

We use additive models to study the yearly evolution of tucuxi dolphin counts, which may be nonlinear, while adjusting for seasonal, water level, and observer effects. Specifically, for each survey we include the year (year), day-of-year (doy), and water level (water) as nonlinear predictors and the number of observers (obs) as a linear predictor. For comparisons, we implement a variety of Bayesian additive models: Gaussian additive models for the raw (am-id) and $\log(y+1)$ transformed data (am-log); star models with identity (am-star-id), unknown parametric (am-star-bc), and unknown nonparametric (am-star-np) transformations; and Poisson (am-Pois) and negative-binomial (am-NB) additive models (using rstanarm). For each method, we jointly evaluate the model performance and the computational efficiency: performance is measured using WAIC, while efficiency is reported as seconds per 10000 effective samples. In particular, we compute multivariate effective sample sizes (Vats et al.,, 2019) for 10 randomly sampled points from the posterior predictive distribution using the mcmcse package in R (Flegal et al.,, 2017). The results are in Table 1. According to WAIC, the star models with unknown transformation (am-star-bc and am-star-np) are strongly preferred, while am-NB is the closest competitor. However, the existing integer-valued additive models (am-Pois and am-NB) are noncompetitive in computational efficiency.

Proceeding with am-star-np, which is selected by WAIC, we plot posterior expectations and credible intervals for each $f_{j}$ in Figure 8. MCMC diagnostics for $f_{j}$ show exceptional mixing with no lack of convergence (see Appendix A.2). The doy plot suggests a seasonal pattern, while the water plot exhibits an approximately quadratic effect. Most importantly, the year plot shows a near linear decline in tuxucis dolphins from 1994-2017, which interestingly has leveled off since 2013. These findings are partially consistent with the results of da Silva et al., (2018) which, assuming only a linear model, also report a significant decrement of dolphins since 1994. Posterior predictive diagnostics are in Appendix A.2, and indicate clear improvements in fit for am-star-np relative to models that exclude transformation (am-star-id) or rounding (am-log).

7 Discussion

star processes provide a mathematically elegant and empirically successful framework to model count and integer-valued data. The approach seamlessly adapts state-of-the-art continuous data models and algorithms to the integer-valued data setting, and thus it offers remarkable modularity both in terms of model specification and computation, while providing ease of implementation and interpretability for practitioners. By incorporating known, unknown parametric, and unknown nonparametric transformations, star processes provide varying degrees of distributional flexibility, and are able to account for important distributional features such as zero-inflation, bounded or censored data, and over- or underdispersion. Empirically, star processes demonstrate goodness-of-fit, out-of-sample point and interval predictive accuracy, reliable inference, and computational scalability.

In addition to the healthcare utilization and animal abundance datasets considered here, we provide further empirical comparisons on three additional datasets in Appendix A.3. Among star, Gaussian, Poisson, and negative-binomial models, the star additive and bart-star models with unknown transformations consistently provide the best performance according to WAIC.

A variety of promising extensions exist for star. The modeling and computational modularity of star suggest that new multivariate, functional, and time series models may be developed for integer-valued data. Furthermore, star is capable of modeling rounded data, which is ubiquitous in practice yet rarely considered in modern statistical and machine learning methods. Lastly, the star model (1)-(2) does not strictly require a Bayesian modeling approach, and may be adapted for classical estimation and inference.

Appendix A Supplemental results

A.1 Evaluating point accuracy for synthetic data

To accompany the WAIC comparisons from the main paper, we evaluate each method for point estimation accuracy. Specifically, we are interested in estimating the conditional expectation of the observed data, $\lambda_{i}^{*}(x)$. For an estimator $\hat{y}_{i}(x)$, we compute the root mean squared error \mbox{RMSE}=\sqrt{\sum_{i=1}^{n}\big{\{}\lambda_{i}^{*}(x)-\hat{y}_{i}(x)\big{\}}^{2}}. The fitted values for star are computed using the conditional expectation of $y$ at $x\in\mathcal{X}$, that is

[TABLE]

where $\mathbb{P}\{y(x)=j\}$ is the star probability mass function and $J(x)$ is a finite truncation. Since $\mathbb{P}\{y(x)=j\}$ depend on the parameters $\theta$ in $\Pi_{\theta}$, the posterior distribution of (14) may be computed by evaluating (14) for each draw of $\theta$ in the MCMC algorithm. Conservatively, we select $J(x)$ to be the 99.99th quantile of the distribution of $y(x)$ pointwise for each $x$, which is easily computable as $h\left[g^{-1}\left\{z_{q}^{*}(x)\right\}\right]$ where $z_{q}^{*}(x)$ is the $q$th quantile of $\Pi_{\theta}$. The point estimate is computed as the posterior expectation of (14).

Figures 9 and 10 depict the relative RMSEs across simulated data sets for the linear and nonlinear simulation designs in Sections 4.1 and 4.2, respectively, defined as the ratio between the RMSE of the generic model over the RMSE for a baseline method, and specifically the lm-log method which represents the common approach of modeling (log-) transformed counts using Gaussian models. Relative RMSE standardizes model performance across simulated datasets: methods with a relative RMSE less than 1.0 demonstrate superior point estimation relative to the baseline method. As in Section 4, we find that star-log and star-bc are consistently competitive and outperform other methods. Interestingly, star-np is much less competitive in RMSE than in WAIC, which suggests that the additional distributional flexibility acquired by modeling $g$ nonparametrically does not necessarily imply more accurate point estimation.

A.2 Model and MCMC diagnostics for the dolphins data

Posterior predictive diagnostics for additive models fit to the tucuxis dolphins data are in Figure 11. The additive star-np model is adequate for the data, while the models which lack either rounding or transformation are incapable of capturing distributional features, including the variability and the proportion of zeros.

The MCMC convergence of the additive star-np model is assessed via traceplots in Figure 12. The traceplots indicate no lack of convergence and demonstrate exceptional mixing: effective sample sizes for all $f_{j}(v_{j})$ exceed 2000.

A.3 Supplemental empirical examples

For further models comparisons, we apply the linear, additive, and bart models to several additional datasets, and again consider star, Gaussian, Poisson, and negative-binomial distributions. For the additive models, each continuous variable (i.e., variables with at least 10 unique observation points) is modeled nonlinearly.

Ships Data

The ships data, available in the MASS package in R, provides the number of damage incidents due to waves for $n=34$ cargo-carrying vessels, as well as ship type (A-E), year of construction (1960-1964, 1965-1969, 1970-1974, or 1975-1979), the period of construction (1960-1974 or 1975-1979), and the aggregated months of service (ranging from 0 to 44882). We model the ship type, year of construction, and period of construction as factors, and center and scale the service variable. The data were analyzed in McCullagh and Nelder, (1989) using a quasi-Poisson regression model to account for observed overdispersion, and subsequently re-analyzed in Mallick and Gelfand, (1994) using a Poisson regression model with unknown link function, which suggests that additional distributional flexibility in the regression model may be important.

Roaches Data

Gelman and Hill, (2006) consider a study of pest management for eliminating cockroaches in city apartments. The response variable, $y_{i}$, is the number of roaches caught in traps in apartment $i$, with $i=1,\ldots,n=262$. A pest management treatment was applied to a subset of 158 apartments, with the remaining 104 apartments receiving a control. Additional data are available on the pre-treatment number of roaches, whether the apartment building is restricted to elderly residents, and the number of days for which the traps were exposed. A notable feature of the data is zero-inflation: $y_{i}=0$ for 94 (36%) of the apartments.

Highway Data

The Highway data, available in the carData package in R, consists of the 1973 accident rate per million vehicle miles on $n=39$ large sections of Minnesota highway. Important predictors include the number of access points per mile, the speed limit, the width of the outer shoulder on the roadway (in feet), and the number of signals per mile of roadway, among others. We consider the accident rate per 10,000 miles, which is the smallest rate for which the observations $y_{i}$ are integer-valued. A notable feature of these data are that, despite being (scaled) accident counts, no two observations $y_{i}$ and $y_{j}$ are equal, and the counts themselves are large, ranging from 161 to 923. Therefore, it is unclear a priori whether an integer-valued model is necessary or advantageous.

Results

The WAICs for the supplementary datasets are reported in Table 2. The star models consistently perform well across all datasets, and in particular star-bc and star-np. Notably, star provides the best linear, additive, and bart models for all datasets with the exception of the Highway data, for which all bart models perform similarly. Interestingly, additive star models are preferred for both the ships data and the Highway data, while bart-star is slightly preferred to the additive star model for the roaches data.