Simultaneous Transformation and Rounding (STAR) Models for Integer-Valued Data
Daniel R. Kowal Antonio Canale


Separation between Second Price Auctions with Personalized Reserves and the Revenue Optimal Auction


    Will MaGoogle Research. willma@google.com.
 Balasubramanian SivanGoogle Research. balusivan@google.com.

    
============================================================================================================






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 () 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 framework, we design a new Bayesian Additive Regression Tree () 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 additive model, which offers greater flexibility and scalability than existing integer-valued models.  The additive model is applied to study the recent decline in Amazon river dolphins.

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




 


















§ 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 <cit.>,  ecology <cit.>,  and insurance <cit.>,  among many others <cit.>.  Counts often serve as an indicator of demand, such as the demand for medical services <cit.>, emergency medical services <cit.>, and call center access <cit.>.
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 <cit.>, negative-binomial <cit.>, zero-inflated Poisson or negative-binomial <cit.>, lognormal Poisson <cit.>, restricted generalized Poisson <cit.>, and Conway-Maxwell Poisson models <cit.>. 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 <cit.>, while log-transformations introduce difficulties in the presence of zeros <cit.>. 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 () 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 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 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 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 models. As a result, provides a cohesive framework for seamlessly adapting state-of-the-art continuous data models and algorithms to the integer-valued data setting.  Using the framework, we design—among others—a new Bayesian Additive Regression Tree () model for integer-valued data. The resulting -model combines the integer-valued distributional flexibility provided by with the predictive and computational advantages inherent to . For synthetic data and a large healthcare utilization dataset (Section <ref>-<ref>), the predictive performance of -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 to study the recent decline in the tucuxi dolphin population, which inhabit the Amazon River.  Using field survey data conducted by <cit.> from 1994 to 2017, we develop a 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 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 <ref> introduces the framework, develops models for the unknown transformation, describes important properties, and discusses computational  details for posterior inference. Section <ref> provides example models, 
which are applied to simulated data in Section <ref> and real data in Sections <ref> and <ref>.

Section <ref> concludes. Additional simulation results and empirical comparisons are in the Appendix. Methods are implemented in the  package  available on GitHub. 







§ SIMULTANEOUSLY TRANSFORMING AND ROUNDING

Consider a count-valued stochastic process y:𝒳→𝒩, where 𝒳 may correspond to predictors, times, or spatial locations and 𝒩 = {0,…,∞}. 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 ℝ or ℝ^+), yet produced a flexible and well-defined distribution on 𝒩. 


To this end, we first introduce continuous-valued process y^*:𝒳→𝒯, 𝒯⊆ℝ related to  the observed count-valued data y via 

    y = h(y^*),

where h: 𝒯→𝒩 is a rounding operator that sets y(x) = j when y^*(x) ∈𝒜_j and {𝒜_j}_j=0^∞ is a known partition of 𝒯. 
For example, we may use the floor function defined by 𝒜_j = [a_j, a_j+1) = [j, j+1) for j ∈𝒩; modifications are available for zero-inflated, bounded, or censored data (Section <ref>).
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 (<ref>). 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 Π_θ on the transformed scale:

    g(y^*) = z^*,    z^*  ∼Π_θ,

where g:𝒯→ℝ is a strictly monotone function. Model (<ref>) is inspired by the common practice of transforming count data prior to application of continuous (Gaussian) models. However, the framework of (<ref>)-(<ref>) defines an integer-valued process for y, in which the transformation g may be modeled as unknown for greater distributional flexibility.









While is sufficiently general to incorporate any continuous family of stochastic process Π_θ for the latent z^*, an important special case of (<ref>) is the conditionally Gaussian regression model:

    z^*(x) = μ(x) + ϵ(x),   ϵ(x) indep∼ N(0, σ^2(x)),

where μ(x) is the conditional expectation of z^*(x) and the errors ϵ(x) are independent but possibly heteroscedastic, conditional on x ∈𝒳. 
Examples of (<ref>) include linear and additive models (Section <ref>) and (Section <ref>), 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  <cit.>.  Our key innovation is the coupling of the transformation (<ref>) with the regression model (<ref>). The transformation g, which we model as unknown, endows the integer-valued process y with greater distributional flexibility, yet leaves model (<ref>) unchanged. This construction allows seamless integration of Bayesian models and algorithms for continuous data of the common form (<ref>)  into the integer-valued framework, with efficient posterior inference available via a general MCMC algorithm (Section <ref>). As demonstrated extensively in the simulations and applications (Sections <ref>-<ref>), models that fail to include both rounding and transformation cannot match the predictive performance of models.










 




The distribution of y is completely determined by the rounding operator h, the transformation g, and the distribution Π_θ. Specifically, the probability mass associated to y(x) = j for each integer j∈𝒩 is 

    ℙ{y(x) = j} = ℙ{y^*(x) ∈𝒜_j} = ℙ{z^*(x) ∈ g(𝒜_j)}.

The distribution of z^* is given by Π_θ, while g(𝒜_j) is determined by the transformation g and the rounding operator h. For model (<ref>) and  𝒜_j = [a_j, a_j+1), 
(<ref>) simplifies to  

    ℙ{y(x) = j} = Φ( g(a_j+1) - μ(x)/σ(x)) -  Φ( g(a_j) - μ(x)/σ(x)).


















The distribution in (<ref>) is related to, yet distinct from, ordinal regression <cit.>. In ordinal regression, each term g(a_j) in (<ref>) is replaced by an unknown latent threshold, say ω_j, with an ordering constraint ω_j ≤ω_j+1 for all j. However, the latent thresholds ω_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 ω_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, 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 <ref>). 




 is fundamentally different from simply rounding the predictions from a continuous data model Π_θ. 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. clearly avoids this issue, and maintains the benefits of using well-known models for continuous data while producing a coherent integer-valued predictive distribution.




 §.§ The transformation g

The transformation g is a crucial component of . When g(t) = t and z^*  is a draw from a Gaussian process,  simplifies to <cit.>. However, the identity transformation is suboptimal in many cases (see Sections <ref> and <ref>). The popularity of log-linear models for count data, especially Poisson and negative-binomial models, suggests that regression effects μ(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) = √(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 <ref>-<ref> 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 θ in the distribution Π_θ of the latent data z^* in (<ref>). Relative to the analogous transformed continuous-valued model, say g(y) ∼Π_θ, the number of parameters is the same, yet produces a coherent integer-valued process.
 A fixed transformation g shares some characteristics with the link function of a generalized linear model (GLM)  <cit.>. For GLMs, the link function maps the expectation of an exponential family distribution to ℝ, which is modeled using a linear predictor. By comparison, maps the continuous-valued y^* to and under (<ref>) models 𝔼[g{y^*(x)} | x] = 𝔼{z^*(x) | x} = μ(x). 


For general application of , 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,  <cit.> similarly relax the assumption of a known link function, adopting a nonparametric approach. For , 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 <cit.>: 

    g(t ; λ) = 
    {(t)|t|^λ - 1}/λ,   λ > 0

with g(t; λ =0) = log(t). Box-Cox functions are a popular choice for transforming continuous data towards Gaussianity, which in the present setting is similar to (<ref>) when Π_θ is Gaussian. Important special cases of (<ref>) include the (shifted) identity transformation g(t; λ = 1) = t-1, the (shifted and scaled) square-root transformation g(t; λ = 1/2) = 2√(|t|) - 2, and the log-transformation. 


To learn the shape of the transformation, we place a prior on λ: we recommend λ∼ 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  <cit.> for g:

    g(t) =  b_I'(t)  γ,

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


We propose a prior for γ in (<ref>) that simultaneously enforces monotonicity, smoothness, and shrinkage toward a pre-specified transformation. However, care must be taken to ensure identifiability of the model and retain interpretability of the parameters θ in Π_θ. For model (<ref>), arbitrary shifting and scaling of g can be matched by shifting and scaling of μ and σ. The parametric transformation (<ref>) preserves identifiability: g(1, λ) = 0  for all λ (shift constraint) and the prior on λ is weakly informative (scale constraint). For nonparametric g in (<ref>), 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→∞ g(t) = 1 (scale constraint), which is enforced by constraining ∑_ℓ=1^L γ_ℓ = 1. Specifically, let

    γ_ℓ = γ̃_ℓ/∑_k = 1^L γ̃_k,   γ̃_ℓindep∼N_+(μ_γ_ℓ, σ_γ^2),   ℓ = 1,…,L,
 
where N_+ is the half-normal distribution. Clearly, γ_ℓ > 0 for each ℓ and ∑_ℓ=1^L γ_ℓ = 1, which guarantees monotonicity of g and preserves identifiability of model (<ref>). We select the prior mean μ_γ = (μ_γ_1,…, μ_γ_L)' such that g(t) is a priori centered around a parametric function of interest, such as (<ref>) with fixed λ = λ_0, and model σ_γ^2  with an inverse-Gamma prior to allow the data to determine the amount of shrinkage toward the parametric function of interest.  Let t_g = (0,1,…, a_max y_i + 1})' and let B_I be the I-spline basis evaluated at t_g, so g( t_g) =  B_I  γ. We solve μ̃_γ = min_μ_γ || g( t_g; λ_0) -   B_I  μ_γ||^2 subject to μ̃_γ_ℓ > 0 for ℓ=1,…,L, which is a one-time cost, 
and normalize μ_γ_ℓ = μ̃_γ_ℓ/∑_k=1^L μ̃_γ_k. In the simulations and applications of Sections <ref>-<ref>, we fix λ_0 = 1/2 and model σ_γ^-2∼(0.001, 0.001).


We use quadratic I-splines with L = 2 + min{( y_i)/4, 10} knots, implemented using the   package in  <cit.>.  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 <cit.> use L=3 or L=5 in all monotone spline examples, a small number of knots may be adequate in many cases. 


















 §.§ Model properties

By design, builds upon models for continuous data, such as those in Section <ref>, and adapts them for integer-valued data. Yet is not merely a mechanism for producing valid integer-valued processes: also provides important distributional properties for modeling integer-valued data in practice. By careful selection of the rounding operator h and the transformation g, 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.  can be parametrized such that zero counts occur whenever 
z^*(x) is negative:

For any model with g(𝒜_0) = (-∞, 0), we have  (i) y(x) = 0 if and only if z^*(x) < 0 and (ii) ℙ{y(x) = 0} = ℙ{z^*(x) ≤ 0}.

Lemma <ref> is valid for known or unknown transformations, and is easily satisfied for (<ref>) letting 𝒜_0 = (a_0, a_1) = (-∞, 1) for λ 0 and 𝒜_0 = (a_0, a_1) = (0, 1) for λ = 0. For model (<ref>), Π_θ is conditionally Gaussian, which may place substantial prior mass on z^*(x) < 0 and thus  y(x) = 0. Therefore, 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: ℙ{y(x) = 0, y(x') = 0} = ℙ{z^*(x) < 0, z^*(x') < 0} depends on the joint distribution of (z^*(x), z^*(x')), which is modeled by Π_θ.





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) ∈{0,1,…,K} and K = 365. can easily incorporate this information into the distribution for y as formalized in the next lemma.

For any model, letting g(𝒜_K) = [g(a_K), ∞) implies  ℙ{y(x) ≤ K} = 1. 

The boundedness constraints in Lemma <ref> are compatible with any choice of (unconstrained) continuous-valued model (<ref>) and do not require modification of the algorithms for estimation and inference in Section <ref>. 
Similar to the case of zero values in Lemma <ref>, allows for dependence among y values that attain the upper bound:  ℙ{y(x) = K, y(x') = K} = ℙ{z^*(x) ≥ g(a_K), z^*(x') ≥ g(a_K)}, which again is modeled by (<ref>) or (<ref>). When K=1, the model (<ref>) with μ(x) = x'β and σ^2(x) = 1 simplifies to probit regression. 


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

Let (h, g, Π_θ) denote model (<ref>)-(<ref>) with rounding operator h, transformation g, and latent distribution Π_θ. 
For right-censored observations y_c(x) = min{y(x), K} with y ∼(h, g, Π_θ) and y_c ∼(h', g, Π_θ) such that h and h' satisfy a_K = a'_K and  𝒜'_K = [a_K, ∞), we have 
ℙ{y(x) ≥ K} = ℙ{z^*(x) ≥ g(a_K)} = ℙ{y_c(x) = K}. 

For right-censored data, the likelihood includes terms of the form ℙ{y(x) ≥ K} for censored observations. Lemma <ref> shows that the censored likelihood terms under a model for y are equivalent to the non-censored likelihood terms ℙ{y_c(x) = K} under a model for y_c = min{y(x), K} with 𝒜'_K = [a_K, ∞). Remarkably, preserves the correct right-censored likelihood for y by directly modeling the observed counts y_c and setting 𝒜_K' = [a_K, ∞), 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 processes are capable of modeling over- or underdispersion. In Figure <ref>, we illustrate the relationships among the expectation 𝔼[y], the variance (y), and the probability of zeros ℙ(y=0) for a process defined by z^* ∼ N(μ,σ^2) and transformation (<ref>) with λ = 1/2. For different values of the parameters the process exhibits different features, including overdispersion, underdispersion, and zero-inflation.  









 §.§ Posterior inference

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

 



Let 𝒟 = {x_i, y_i}_i=1^n denote the observed pairs of points x_i ∈𝒳 and integer-valued data y_i = y(x_i). Consider a Bayesian specification of (<ref>) with suitable prior on  θ  and an algorithm 𝔸 which draws from the posterior distribution of θ given the (continuous) data.  The sampling algorithm 𝔸 is designed for continuous data, such as (Gaussian) additive models or , depending on the different choices for Π_θ (see Section <ref> for details). 
The posterior sampling algorithm for defines a Gibbs sampler by combining a data augmentation step with algorithm 𝔸 as follows:


  * Sample [z^*(x_i) |𝒟, θ] from Π_θ truncated to  z^*(x_i) ∈ g(𝒜_y_i) for i =1, … n; 

  * Sample [θ| z^* ] using algorithm 𝔸 conditioning on z^* = (z^*(x_1),…, z^*(x_n))'.

In the case of model (<ref>), the data augmentation step may be computed efficiently using a standard univariate truncated normal sampler. 
Specifically, for 𝒜_j = [a_j, a_j+1), the full conditional distribution of the latent data is [z^*(x_i) |𝒟, θ] ∼ N(μ(x_i), σ^2(x_i)) truncated to [g(a_y_i), g(a_y_i + 1)). 
While the process y^* is useful for interpretability of the 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 (<ref>) this translates to sampling the parameter λ from its full conditional posterior distribution, for which we use a slice sampler <cit.>. For the nonparametric model (<ref>) with prior (<ref>), we sample ξ_γ = log(γ̃) using Metropolis-Hastings and [σ_γ^-2| - ] ∼{0.001 + L/2, 0.001 + ∑_ℓ=1^L (γ̃_ℓ - μ_γ_ℓ)^2}, and set g(t) = b_I'(t)γ as defined in (<ref>)-(<ref>). 

The sampler for ξ_γ uses a Gaussian random walk proposal with covariance matrix tuned using the robust adaptive Metropolis (RAM) algorithm of <cit.> 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  <cit.> 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 [ỹ(x) | 𝒟], we additionally sample [z̃^*(x) |θ] from Π_θ using the current draw of θ and set ỹ(x) =h[g^-1{z̃^*(x)}] for each x. This step is extremely simple, yet provides inference for integer-valued predictions, model-based imputation of missing data at  x ∈𝒳, and useful model diagnostics.  For parametric g  in (<ref>), the functions g^-1(s; λ) are known, while for nonparametric g we approximate g^-1(s)  ≈min_t | s  -  b_I'(t) γ|, 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 <cit.> 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 <cit.>, factor models <cit.>, and functional time series models <cit.>, yet faces two important limitations: first, it is restricted to the negative-binomial distribution, and second, the resulting MCMC sampler is often inefficient  <cit.>. As demonstrated in Section <ref> and Appendix <ref>, the proposed algorithm provides excellent MCMC efficiency, even for nonlinear versions of model  (<ref>). 




§ REGRESSION MODELING WITH STAR

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















 §.§ Linear regression

The simplest and most common special case of (<ref>) is linear regression, where σ^2(x) = σ^2 and μ(x) = x'β for x a p-dimensional vector of predictors. Let X denote the n × p matrix of predictors, and consider the (independent) priors [β] ∼ N(0, Σ_β) and [σ^-2] ∼(a_σ, b_σ), where Σ_β is the prior covariance for the regression coefficients. We focus on the ridge prior [β_j | σ_β] iid∼N(0, σ_β^2) for illustration, with [σ_β] ∼(0, A), but extensions for other shrinkage priors are available <cit.>. The Gibbs sampler for the linear model is given by the following steps:
 

  * Sample [z^*(x_i) | -] ∼ N(x_i'β, σ^2) truncated to  z^*(x_i) ∈ g(𝒜_y_i) for i =1, … n; 



  * Sample [β| -] ∼ N(Q_β^-1ℓ_β, Q_β^-1) where Q_β = σ^-2 X'X +   Σ_β^-1 and ℓ_β = σ^-2 X' z^*;


  * Sample [σ^-2|  - ] ∼[a_σ + n/2, b_σ + ∑_i=1^n {z^*(x_i) - x_i'β}^2/2];


  * Sample [σ_β^-2|  - ] ∼{(p-1)/2, ∑_j=1^p β_j^2} truncated to [1/A^2, 0) and set Σ_β = σ_β^2  I_p.

 
For other conditionally Gaussian shrinkage priors, only the sampler for [Σ_β | -] in the final step must be modified. Note that we include a diffuse prior on the intercept, β_0 ∼ N(0, 10^6), and select the hyperparameters a_σ = b_σ = 0.001 and A = 10^4 for Sections <ref>-<ref>.






Computational scalability of the linear model is primarily determined by the sampler for the regression coefficients β, which depends on both the sample size n and the number of predictors p. Importantly, computationally efficient samplers for [β| z^*, σ^2] exist for p<n <cit.> and p≫ n <cit.> for conditionally Gaussian priors. 
By design, the MCMC algorithm leverages these existing priors and algorithms for Gaussian data to provide efficient and scalable posterior inference for integer-valued data. 















 §.§ Additive models

Suppose the p predictors are partitioned as x' = (u', v') for linear predictors u and nonlinear predictors v. The additive model is given by (<ref>) with conditional mean

    μ(x) = u'β + ∑_j f_j(v_j),

where f_j :𝒳_j→ℝ is an unknown function of v_j ∈𝒳_j. The unknown f_j are typically modeled as smooth nonparametric functions, and may capture nonlinearities in each v_j. The linear model is a special case of (<ref>) with μ(x) = x'β. For the conditional variance of z^* in (<ref>), we use the conditionally conjugate prior σ^-2∼(0.001, 0.001). 


Within the framework, we apply flexible and computationally efficient parametrizations of (<ref>) that have been well-developed for Gaussian and  exponential family models. The linear regression coefficients are assigned conditionally Gaussian priors, β∼ N(0, Σ_β), including the ridge prior Σ_β = σ_β^2 I and other shrinkage priors <cit.> as special cases. The nonlinear functions in (<ref>) are modeled smoothly using a basis expansion f_j(v_j) =  b_j'(v_j)  α_j, where b_j is an L_j-dimensional vector of known basis functions and α_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 α_j ∼ N(0, σ_α_j^2 I) is the smoothing prior and B_j'B_j is diagonal, where B_j = (b_j(v_1,j), …, b_j(v_n,j))' is the basis matrix <cit.>. The nonlinear terms are constrained such that ∑_i=1^n f_j(v_i,j) = 0 for identifiability, which is enforced in the reparametrization of the basis b_j. We let σ_α_j^-2∼(0.1, 0.1), which allows the smoothness of each f_j to be learned from the data. 


For observed predictors x_i' = (u_i', v_i'), let U denote the matrix of linear predictors and let f_j =  B_j  α_j. The Gibbs sampler for the additive model iterates the following full conditional distributions:
 


  * Sample [z^*(x_i) | -] ∼ N(u_i'β + ∑_j f_j(v_i,j), σ^2) truncated to  z^*(x_i) ∈ g(𝒜_y_i); 



  * Sample [β| -] ∼ N(Q_β^-1ℓ_β, Q_β^-1) where Q_β = σ^-2 U'U +   Σ_β^-1 and ℓ_β = σ^-2 U'( z^* - ∑_j  f_j);



  * For each j, sample [ α_j | -] ∼ N(Q_α_j^-1ℓ_α_j, Q_α_j^-1) where Q_α_j = σ^-2 B_j'B_j +   σ_α_j^-2I and ℓ_α_j = σ^-2 B_j'( z^* - Uβ - 
∑_k  j f_k)  and set f_j = B_j α_j;




  * Sample [σ^-2|  - ] ∼(0.001 + n/2, 
0.001 + || z^*  -  U β - ∑_j  f_j||^2/2);




  * For each j, sample [σ_α_j^-2| -] ∼(0.1 + L_j/2, 0.1 + ∑_ℓ=1^L_jα_j^2/2).

 
The MCMC sampling algorithm for the additive model is efficient, with empirical support provided in Section <ref>. The computational complexity for the nonlinear basis coefficients {α_j} is 𝒪(∑_j L_j) due to the diagonality of B_j'B_j, while the linear coefficients β also may be sampled efficiently <cit.>.  Additional sampling steps for Σ_β depend on the model specification, but are often available in closed form. 

























 §.§ 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, <cit.> proposed , which is a “sum-of-trees" model within a fully Bayesian framework. Tree-based regression models, such as <cit.>, are designed to model complex interactions among predictors. Notably, utilizes many trees, where each individual tree is constrained via the prior to be a weak learner. As a result, 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 is highly competitive with state-of-the-art statistical and machine learning models. 

For integer-valued data, has been relatively underutilized. Adaptations of for negative-binomial data are feasible via Pólya-Gamma augmentation <cit.>, similar to the probit implementation in <cit.> 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) sampling algorithm. Recently, <cit.> proposed a log-linear 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 framework, we parametrize the model (-) as in 
<cit.> and specifically

    μ(x) = ∑_k=1^m f(x; T_k, M_k),

where T_k is a binary tree comprised of interior splitting rules and terminal nodes and  M_k = {η_1,k,…, η_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 (<ref>) assigns a value η_ℓ,k∈ M_k, and these values are summed across all trees k=1,…,m. <cit.> propose prior distributions that constrain each T_k to be shallow, thereby limiting the order of interactions, and constrain each η_ℓ,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 <cit.>, and a prior for the terminal values given the tree, p(η_ℓ,k | T_k), which is Gaussian.  A key feature of  is that, by transforming to Gaussianity, we inherit the same framework as the original , and therefore may directly incorporate the well-studied priors and hyperparameters from <cit.>. 


More careful consideration is required for the prior distribution of  σ^2 in (<ref>).  <cit.> emphasize that an informative prior distribution is important to balance between overly aggressive and overly conservative model fits, and parametrize the prior for σ^2 as an inverse chi-square distribution calibrated using a data-based overestimate  σ̂  of σ. However, any statistics calculated from the data y are likely inappropriate for , since the transformation g impacts the scale of z^*. As a remedy, we compute σ̂ as the posterior median of σ from the linear model (<ref>), where the transformation in the linear model is chosen to match the transformation in -. Given σ̂, which indeed is a data-based overestimate of σ, we adopt the default hyperparameter suggestions of <cit.>. 


For posterior inference under -, we combine a data augmentation step for z^*(x_i) similar to  Section <ref> with a sweep from the original <cit.> algorithm to draw the parameters in (<ref>) using z^* as data. The <cit.> sampler proceeds using backfitting: draws for the kth tree [(T_k, M_k) |   z^*, {(T_k', M_k')}_k  k', σ] are generated using <cit.> and σ^2 is sampled from an inverse-Gamma distribution. Incorporating these sampling steps  into the larger Gibbs sampler in Section <ref> is straightforward using the  package <cit.>. 






§ SIMULATION STUDIES

The proposed 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, … n and n=100 are simulated from a negative-binomial distribution parametrized by conditional mean 𝔼{y_i(x) | λ_i^*} = λ_i^*(x) and variance {y_i(x) | λ_i^*, r^*} = λ_i^*(x){1 + λ_i^*(x)/r^*} with dispersion parameter r^* > 0. 




As r^* decreases to zero, the variance increasingly dominates the mean while as r^*→∞, the distribution converges to a Poisson distribution with parameter λ_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 <ref>) and nonlinear (Section <ref>) parametrizations for the log-mean logλ_i^*(x).

We emphasize that in all cases, the simulated datasets are not generated under the proposed model: they are simulated from negative-binomial and (approximate) Poisson distributions.


Competing models are compared using the Watanabe-Akaike/widely-applicable information criteria (WAIC)  <cit.>. 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 ℳ is defined as _ℳ = -2(_ℳ - d_ℳ), 
where d_ℳ is the effective number of parameters for model ℳ and _ℳ is the log-predictive pointwise density defined by 

    _ℳ(y) = ∑_i=1^n log( 1/S∑_s=1^S p_ℳ(y_i |θ^s) )


for θ^s drawn from its posterior distribution. For models, we simply have


    _(y) = ∑_i=1^n log( 1/S∑_s=1^S  
    Φ{g^s(a_y_i+1) - μ^s(x_i)/σ^s(x_i)} - 
    Φ{g^s(a_y_i) - μ^s(x_i)/σ^s(x_i)}).

For the effective number of parameters, we follow the recommendation of <cit.> and use the sample variance of the pointwise log-likelihoods across MCMC simulations: 

d_ℳ  = ∑_i=1^n ( log p(y_i |θ^s)).

The pointwise log-likelihood of is simple and efficient to compute, and is sufficient for computing WAIC as well as other information criteria.



WAIC is used in Sections <ref> and <ref> for synthetic data, but exact out-of-sample metrics are provided for the real data in Section <ref>. Appendix <ref> also provides root mean square errors for estimating the conditional expectation λ_i^*(x) to evaluate point estimation accuracy. 





 §.§ Linear mean functions

We first consider a linear log-mean,  logλ_i^*(x) = β_0 + ∑_j=1^p x_i,jβ_j, where the p=6 predictors are drawn independently from x_i,j∼ N(0,1) and the coefficients are  β_0 = log(1.5), β_1 = β_2 =β_3 = log(2.0), and   β_4=β_5=β_6 = 0. Under this specification, the expected counts at x_i,j = 0 is 1.5, while each nonzero coefficient β_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 models, we use the linear model (<ref>) with prior [β_j | σ_β] ∼ N(0, σ_β^2) and [σ_β] ∼(0, 10^4), and the following transformations: known transformation (<ref>) with λ = 0 (lm-star-log), λ = 1/2 (lm-star-sqrt), and λ = 1 (lm-star-id); unknown parametric transformation (<ref>) (lm-star-bc); and unknown nonparametric transformation (<ref>) (lm-star-np). We also include the same linear model, but  instead with the Gaussian model (<ref>) applied directly to the counts y (lm) and the log-transformed counts log(y + 1)  (lm-log). These models are natural competitors to , since they incorporate the same model for μ(x) but omit the rounding step in (<ref>) 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  package <cit.>. lm-Pois and lm-NegBin are widely used for modeling count data, and here correspond to the true data-generating process. 

In Figure <ref>, 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 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 model improvements relative to  the Gaussian models and the identity transformation lm-star-id definitively demonstrate the importance of both rounding and transformation.















 §.§ Nonlinear mean functions

To evaluate bart-star, we specify a nonlinear form for the log-mean,   logλ_i^*(x) = β_0 + β_1 f̃(x), where  f̃(x) is the centered and scaled Friedman function <cit.>

    f(x) = 10 sin(π x_1 x_2) + 20(x_3 - 0.5)^2 + 10x_4 + 5x_5

featured in the original simulations  <cit.>. As in <cit.>, we select p=10 and simulate x_i,jiid∼(0,1). We fix the parameters β_0 = log(1.5) and β_1 =log(5.0), which again corresponds to low counts with a moderate signal. 


We combine the -model of Section <ref> with known transformation (<ref>) for λ = 0 (bart-star-log), λ = 1/2 (bart-star-sqrt), and λ = 1 (bart-star-id); unknown parametric transformation (<ref>) (bart-star-bc); and unknown nonparametric transformation (<ref>) (bart-star-np). For competitors, we include the Gaussian model (-id) of <cit.> and a Gaussian model on the log-transformed counts log(y+1) (-log). Lastly, we include the linear models lm-star-bc and lm-log from Section <ref>. 


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

















§ 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 models for this task, we use data from the National Medical Expenditure Survey (NMES)  available in the  package in  <cit.>. Multiple measures of healthcare utilization are available, including physician office visits (), non-physician office visits (), physician hospital outpatient visits (), and non-physician hospital outpatient visits (). 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  <cit.> and  <cit.>. 






The NMES data provides a unique opportunity for insightful out-of-sample prediction comparisons. Each measure of healthcare utilization  (, , , and ) is count-valued with distinct characteristics: the probability mass functions in Figure <ref> 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 <ref>. 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 <ref>; for conciseness, we report results for -log, --id, --log, and --bc. The (untransformed) model of <cit.> 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 {ỹ_i}_i=1^n_test for each model ℳ:

    _ℳ(ỹ) = 1/n_test∑_i=1^n_testlog p_ℳ(ỹ_i | y) ≈1/n_test∑_i=1^n_testlog{1/S∑_s=1^S p_ℳ(ỹ_i |θ^s) }

where θ^s ∼ p_ℳ(θ | y) is a draw from the posterior under model ℳ. For with model (<ref>), we have 
p_ℳ(ỹ_i |θ^s) = Φ[ {g^s(a_ỹ_i  +1) - μ^s(x_i)]}/σ^s(x_i)] -  Φ[ {g^s(a_ỹ_i) - μ^s(x_i)}/σ^s(x_i)] similar to (<ref>). 

The results for each response variable are in Figure <ref>; larger values indicate more accurate predictive distributions. The model with unknown transformation (--bc) performs best in all cases, closely followed by the model with known transformation (--log). Omitting the transformation (--id), similar to <cit.>, leads to substantial deterioration in predictive accuracy. Furthermore, the untransformed --id model produced infinite scores for some test points ỹ_i; these points were excluded for computing (<ref>), but imply that --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 {ỹ_i}_i=1^n_test: smaller intervals that provide the correct coverage are preferable. The MPIWs with empirical coverages are displayed in Figure <ref>. Across all responses, the model with unknown transformation (--bc) provides the most precise prediction intervals with  correct coverage. For the responses with greater zero-inflation and overdispersion (, , and ), --bc reduces MPIWs by a median of 45%, 50%, and 59%, respectively, relative to --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, ℙ(ỹ_i > 0 | y). We evaluate each method using logarithmic scoring on the event {ỹ_i > 0}, which is a proper scoring rule for binary outcomes <cit.>. The results are in Figure <ref>, where larger values indicate superior predictive accuracy. The -models which include both transformation and rounding (--log and --bc) decisively outperform the competitors, especially the continuous-valued model (-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 models, and in particular -models which include both transformation and rounding. 




§ 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, <cit.> 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 <cit.> fit a linear model to the logarithm of dolphin counts, we propose to leverage the 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 (), day-of-year (), and water level () as nonlinear predictors and the number of observers () 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); models with identity (am--id), unknown parametric (am--bc), and unknown nonparametric (am--np) transformations; and Poisson (am-Pois) and negative-binomial (am-NB) additive models (using ). 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 <cit.> for 10 randomly sampled points from the posterior predictive distribution using the  package in  <cit.>.  The results are in Table <ref>. According to WAIC, the models with unknown transformation (am--bc and am--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--np, which is selected by WAIC, we plot posterior expectations and credible intervals for each f_j in Figure <ref>. MCMC diagnostics for f_j show exceptional  mixing with no lack of convergence (see Appendix <ref>). 
The  plot suggests a seasonal pattern, while the  plot exhibits an approximately quadratic effect. Most importantly, the  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 <cit.> which, assuming only a linear model, also report a significant decrement of dolphins since 1994. Posterior predictive diagnostics are in Appendix <ref>, and indicate clear improvements in fit for am--np relative to models that exclude transformation (am--id) or rounding (am-log). 












§ DISCUSSION

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, 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, 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 <ref>. Among , Gaussian, Poisson, and negative-binomial models, the  additive and -models with unknown transformations consistently provide the best performance according to WAIC.

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







§ SUPPLEMENTAL RESULTS
 



 §.§ 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, λ_i^*(x). For an estimator ŷ_i(x), we compute the root mean squared error = √(∑_i=1^n {λ_i^*(x) - ŷ_i(x)}^2). The fitted values for are computed using the conditional expectation of y at x ∈𝒳, that is

    𝔼{y(x)} = ∑_j=0^∞ j ℙ{y(x) = j}≈∑_j=1^J(x) j ℙ{y(x) = j},

where ℙ{y(x) = j} is the probability mass function and J(x) is a finite truncation. Since ℙ{y(x) = j} depend on the parameters θ in Π_θ, the posterior distribution of (<ref>)
may be computed by evaluating (<ref>) for each draw of θ 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[g^-1{z_q^*(x)}] where z_q^*(x) is the qth quantile of Π_θ. The point estimate is computed as the posterior expectation of (<ref>). 

Figures <ref> and <ref> 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 -log and -bc are consistently competitive and outperform other methods. Interestingly, -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. 





















 §.§ Model and MCMC diagnostics for the dolphins data














Posterior predictive diagnostics for additive models fit to the tucuxis dolphins data are in Figure <ref>. The additive -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 -np model is assessed via traceplots in Figure <ref>. The traceplots indicate no lack of convergence and demonstrate exceptional mixing: effective sample sizes for all f_j(v_j) exceed 2000. 














 §.§ Supplemental empirical examples

For further models comparisons, we apply the linear, additive, and models to several additional datasets, and again consider , 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  data, available in the  package in , 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 <cit.> using a quasi-Poisson regression model to account for observed overdispersion, and subsequently re-analyzed in <cit.> using a Poisson regression model with unknown link function, which suggests that additional distributional flexibility in the regression model may be important. 

















 











  §.§.§ Roaches Data

<cit.> 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,…, 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  data, available in the  package in ,  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 <ref>. The models consistently perform well across all datasets, and in particular -bc and -np. Notably, provides the best linear, additive, and models for all datasets with the exception of the  data, for which all models perform similarly. Interestingly, additive models are preferred for both the  data and the  data, while -is slightly preferred to the additive model for the  data.
















 
 
  
  
 
 
 






apalike




