Suppose you have a Bayesian statistical model for continuous responses that you believe works extremely well in several settings. Refer to this statistical model as the "preferred model." Also suppose you have observed a data set consisting of multiple response types (e.g., continuous, countvalued, Bernoulli trials etc.). These response types may be "mismatched" with the response types of the preferred model. For example, the data set may consist of count-valued observations, but this preferred model may be derived only for Gaussian data. The primary goal of this article is to introduce a reasonable easy-to-implement all-purpose method that "converts" a Bayesian statistical model for continuous responses into a Bayesian model appropriate for the analysis of multi-response data sets.
There are several methods for jointly modeling data consisting of multiple response types, however these approaches require one to either abandon the preferred model, or it requires you modify it in a manner that creates computational difficulties. For example, multi-task learning models (Argyriou et al., 2007; Kim and Xing, 2009; Yang et al., 2009) , regression trees, and random forests (Hastie et al., 2009; Fellinghauer et al., 2013) can jointly model multiple response types, but does not allow an analyst to incorporate their preferred model for continuous data. Yang et al. (2014) use a specific model (that may differ from the preferred model) and use a non-Bayesian approach for implementation. Copulas (Liu et al., 2009; Xue and Zou, 2012; Dobra and Lenkoski, 2011; Liu et al., 2012) are also available for analyzing such data, but again, they do not immediately incorporate an analyst's preferred model.
A standard approach to model non-Gaussian data (often of a single type) with a preferred model is to define a generalized linear mixed effects model (GLMM) (e.g., see McCulloch et al., 2008 , for a standard reference). For example, Bernoulli data is modeled hierarchically, where the logit of the probability of success can be analyzed using an analysts preferred model. GLMMs lack conjugacy, which creates noticeable difficulty when implementing a GLMM on a modern high-dimensional data set. A more recent alternative is the LCM. Basic theoretical results and empirical analyses in Bradley et al. (2018) , Hu and Bradley (2018) , H.-C. Yang et al. (2019) , Bradley et al. (2019c) , and Bradley et al. (2019a) suggest that one can outperform a standard GLMM (specifically Latent Gaussian Process (LGP) models) in terms of prediction error. However, the LCM requires the preferred model to be a mixed effects model and to modify the distribution of random effects to follow the appropriate distribution based on conjugacy.
A classical approach that is less commonly used (in a formal context) for Bayesian models is to transform the data, so that the transformed data can be reasonably modeled using the distribution assumed by the preferred model. In the non-Bayesian setting this literature is extremely well-developed and includes the Box-Cox transformations (Box and Cox, 1964) , the alternating conditional expectations (ACE; Breiman and Friedman, 1985) algorithm, graphical techniques (McCulloch, 1993) , and the Yeo-Johnson power transformation (Yeo and Johnson, 2000) , among other techniques. More recently developments in rank based algorithms (Servin and Stephens, 2007; McCaw et al., 2019; Beasley et al., 2009 ) and quantile-matching (McCullagh and Tresoldi, 2020) have also been proposed in the non-Bayesian setting. It is important to note that Bayesian models for transformations have been proposed as well, but focus on the case where continuous non-normal data are observed and the preferred model assumes normality. In particular, these Bayesian models put a prior on the free parameter within the Box-Cox transformation or the Yeo-Johnson power transformation (Kim et al., 2013; Charitidou et al., 2015 Charitidou et al., , 2018 . No such Bayesian model has been developed to analyze multi-response response data using any preferred model for a continuous response.
There are three distributions that define our proposed model: (a) the distribution of the data given a transformation, (b) the prior distribution of the transformation, and (c) the distribution of the process of interest given the transformation (i.e., the aforementioned preferred model). In this article, we model the data given a transformation (a) using members from the exponential family.
Specifically, given a transformation, continuous data follows the normal distribution, categorical data follows the multinomial distribution, binary data follow the Bernoulli distribution, and countdata follow the Poisson distribution. These distributions are conjugate with the logit-beta (Gao and Bradley, 2019; Bradley et al., 2019d) and the log-gamma distributions (Bradley et al., 2018; Hu and Bradley, 2018; Bradley et al., 2019a; H.-C.Yang et al., 2019) , which are special cases of the Diaconis-Ylvishaker (DY) distribution (e.g., see Diaconis and Ylvisaker, 1979; Chen and Ibrahim, 2003 , for key references). Consequently, the prior distribution of the transformation (b) is modeled with a DY distribution. This is a new way in Bayesian statistics to model non-Gaussian processes using models for continuous data. The implementation of our approach is straightforward, and can be done using composite sampling. In particular, the first step is to sample from the posterior distribution of the transformation. Then the second step is to sample from the conditional distribution of the latent process of interest given the transformation. This conditional distribution is derived from the preferred model. The first step of the composite sampler is computationally straightforward because the DY distribution is conjugate (and easy to sample from) with the exponential family. Additionally, the first step of this algorithm is important for the purpose of analyzing multiple response types.
Specifically, at the end of the first step we obtain a replicate from the posterior distribution of the transformation (which is continuous valued). Thus, the first step of the composite sampling algorithm "transforms" the multi-response data into a continuous-valued quantity appropriate for the preferred model.
Implementation of the preferred model is unchanged in the second step of our composite sampling algorithm. This is particularly noteworthy, as many of the Bayesian statistical models derived for Gaussian data are not immediately computationally efficient in the non-Gaussian data setting (e.g., see Bradley et al., 2019b; Kang and Cressie, 2011; Katzfuss and Cressie, 2012, for examples in the spatial setting). This is because GLMMs in the non-Gaussian setting have full-conditional distributions that are not Gaussian, and can not be sampled from immediately. Bayesian meth-ods that do not have easy to sample from full-conditional distributions require difficult to tune Metropolis-Hastings algorithms (e.g., see Bradley et al., 2019a , for an example), inefficient rejection samplers (e.g., see Damien et al., 1999) , or significant reparameterization to make approximate Bayesian methods (that are only appropriate for small parameter spaces) practical (Rue et al., 2009; Neal, 2011) . The second step of our composite sampling algorithm allows one to circumvent this issue entirely, and simply use the computational strategies that were developed for the preferred model.
The two steps of our composite sampler can be seen as sequential smoothing. By "smoothing"
we mean a function of the data that attempts to discover important features in the data (e.g., see Simonoff, 2012 , for a standard reference). Multiple layers of smoothing may lead to estimates that are "oversmooth," in the sense that many features of the data are not captured. To avoid oversmoothing we specify the model so that the posterior distribution of the transformation is "saturated." Recall a saturated model is one in which there exists at least as many parameters as there are data points, and fitting this model allows you to exactly recover the original data set.
Hence, saturated models are often an extreme example of overfitting. Thus, in the first step of our composite sampler we choose to overfit the data, and in the second step we smooth overfitted values (again this is done to avoid oversmoothing).
In the classical log-linear model literature, saturated models are useful for selecting more parsimonious models (e.g., see Agresti, 2007 , for a standard reference). Specifically, the most parsimonious reduced model that is not significantly different (in terms of the deviance or chi-square statistic) from the saturated model is used for statistical inference. Consequently, specifying the transformation model to be saturated allows us to assess the goodness-of-fit of the preferred model in a fully Bayesian manner that is similar to what is done in classical residual analysis.
The remainder of this article is organized as follows. In Section 2, we review transformations of data. Then, we introduce general Bayesian methodology to analyze multi-response data with unknown transformations. Additionally, we provide a specific class of transformation models and describe goodness of fit. In Section 3, we give simulation studies to show that our approach can be used in several methodological settings. In particular, we apply our approach to BART models and SME models. Additionally, in Section 4 we provide a joint spatial analysis of a spatiotemporal zero-inflated data set (i.e., continuous and binary responses) consisting of public costs of natural disasters (in billions of USD) records obtained from the National Oceanic Atmospheric Association's (NOAA) National Centers for Environmental Information (NCEI). To analyze this data we apply our approach to the multivariate spatio-temporal mixed effects model (MSTM) from Bradley et al. (2015) , which is a model that allows for non-stationary temporal dynamics and flexible multivariate spatial dependence for Gaussian data.
Let Z i j denote the observed data for i = 1, . . . , I j and j = 1, 2, 3. We consider the setting where for each i, Z i1 is continuous-valued, Z i2 is integer-valued ranging from 0, . . . , b i , and Z i3 is binary. It is assumed that Z i1 , Z i2 , and Z i3 are distributed according to different classes of probability density functions/probability mass functions (hence referred to as multiple-type responses). One classical strategy to model data of this type is to impose a transformation,
where h(·) is a transformation of the datum Z i j , "Dist" is a short-hand used for a probability density
is known as a link function (e.g., see
McCullagh and Nelder, 1989, for a standard reference). Additionally, Y i j is defined for i = 1, . . ., I and j = 1, 2, 3, where I ≥ max(I 1 , I 2 , I 3 ). Here, "Dist Y i j , θ θ θ " represents the aforementioned preferred model. In what remains, inference on {Y i j } and θ θ θ is the primary goal.
To aid in our exposition we drop the functional notation for h(·) and write h i j = h j (Z i j ). Transformations convert a multiple response type data set (e.g., {Z i j }) to a single response type data set (e.g., {h i j }), since h i j follows a single distribution with a continuous support. Consequently, transformations have become a standard tool in analyzing multiple response types. Recall, transformations such as these have a long history including the box-cox transformations (Box and Cox, 1964) , graphical techniques (McCulloch, 1993) , the alternating conditional expectations (ACE; Breiman and Friedman, 1985) algorithm, and the Yeo-Johnson power transformation (Yeo and Johnson, 2000, among others).
In this paper, we introduce a Bayesian solution to the problem of an unknown transformation.
In particular, we define pdf and probability mass functions (pmf), f (Z i j |h i j ). We refer to these distributions as "transformation models." In Section 2.2, we describe Bayesian implementation (Section 2.2) using a general transformation model and any well defined preferred model. Then, in Section 2.3 an example specification of the transformation model is given. Finally we discuss goodness-of-fit in this Bayesian paradigm in Section 2.4.
In this section, we describe Bayesian implementation of the model introduced in Section 2.1.
Here, let n = ∑ 3 j=1 I j , the n-dimensional data vector z = (Z 11 , . . . , Z I 3 3 ) ′ , the n-dimensional trans-
From (1), the preferred model "Dist" is represented in terms of a hierarchical model:
Following the terminology used in Cressie and Wikle (2011), we call f (h i j |Y i j , θ θ θ ) the "transformed data model," f (Y i j |θ θ θ ) the "process model," and f (θ θ θ ) the "parameter model" (or prior). Bayes rule can be used to produce the following conditional distribution (e.g., see Gelman et al., 2013 , for a standard reference),
where we have assumed h i j is conditionally independent of h km given Y i j and θ θ θ for k = i and m = j
Similarly, one can use Bayes rule to produce the posterior distribution of the transformed data. That is,
where the distribution f (h|γ γ γ) is referred to as a "transformation prior," the q-dimensional realvalued vector γ γ γ is referred to as a transformation hyperparameter, and the distribution f (γ γ γ) is referred to as a "transformation hyperprior."
Equations (2) and (3) can be used to produce a posterior distribution for y and θ θ θ .
is conditionally independent of γ γ γ given h, and z and (y ′ , θ θ θ ′ ) ′ are conditionally independent given h. Then:
Proof: See Appendix A.
The model in (4) can easily be simulated from using a composite sampling scheme, provided that it is easy to simulate from f (y, θ θ θ |h). Algorithm 1 gives the step-by-step implementation of how to Algorithm 1: Steps Needed for Simulating from (4).
1: Set b = 1 and initialize h, γ γ γ, y, and θ θ θ with h [0] , γ γ γ [0] , y [0] , and θ θ θ [0] .
from their full-conditional distributions. We use the slice sampler (Neal et al., 2003) if the full-conditional does not have a closed form. 4: Sample y [b] and θ θ θ [b] from f (y, θ θ θ |h [b] ), which is the posterior distribution associated with the preferred model described in (2).
simulate from the posterior distribution in (4). Here, we see that the implementation of our model is similar to the bootstrap implementation, where we have replaced a resampling step with sampling from f (h|z) and the full-conditional distributions associated with γ γ γ. This similarity emphasizes the flexibility of allowing for unknown transformations in a Bayesian context, since the bootstrap algorithm is an established flexible approach in the literature (e.g., see Efron, 1992 , for an early reference). Of course, the bootstrap algorithm produces replicates from a different distribution than that of Algorithm 1. Specifically, the bootstrap method results in an approximate sample from the sampling distribution of a statistic. Whereas, the composite sampling approach in Algorithm 1 can be seen as a means to sample from (4). This is also different from the Bayesian bootstrap (Rubin, 1981) , which does not restrict the samples to be from a posterior distribution of the form in (4).
Consider the following specifications of the transformation model:
where Normal(h i1 , v) is a shorthand for the normal distribution with mean h i j ∈ R and variance v > 0, Binomial(b i , p) is a shorthand for the binomial distribution with b i > 1 number of trials and probability of success p ∈ (0, 1), and Poisson(µ i j ) is a shorthand for the Poisson distribution with
With the transformation model f (z|h) defined, we are left to specify a transformation prior and transformation hyperprior. We define the transformation prior to be the conjugate distributions associated with (5). It follows from Diaconis and Ylvisaker (1979) that the conjugate distribution for h i j is given by,
where
log(x). Also, we use the shorthand DY(α j , κ j ; ψ j ) to represent the pdf in (6). Finally, let γ γ γ = (α 1 , α 2 , α 3 , κ 1 , κ 2 , κ 3 , v) ′ be the transformation hyperparameter.
Equations (5) and (6) can be used to produce a full-conditional distribution for h.
Proposition 2: Suppose (5) and (6) holds. Then,
Proof: See Appendix A.
It is straightforward to simulate directly from the posterior distribution in (7). Replicates of h i j from (7) can be computed using the following transformation (Bradley et al., 2019a) :
where " d =" stands for equal in distribution, w 1 |Z i1 , α 1 , κ 1 , v is distributed normally with mean zero and variance 2κ 1 + 1 v −1 , w 2 |Z i2 , α 2 , κ 2 is distributed according to a beta distribution with shape parameters (α 2 + Z i2 ) and (κ 2 − α 2 + b i − Z i2 ), and w 3 |Z i3 , α 3 , κ 3 is distributed according to a gamma distribution with shape parameter (α 3 + Z i3 ) and rate parameter (κ 3 + 1).
Step 2 of Algorithm 1 involves simulating according to (8).
The specification of a transformation hyperprior for γ γ γ is crucial to guarantee that f (h i j |Z i j , γ γ γ) is
proper in the event that Z i3 = 0, Z i2 = 0, or Z i3 = b i . Thus, we assume α 1 = κ 1 = 0, α 2 and α 3 are distributed according to a gamma distribution, κ 2 |α 2 is distributed according to a shifted (by α i ) gamma distribution, κ 3 follows a gamma distribution, and v is distributed according to an inverse gamma distribution (e.g., see Gelman, 2006, among others). These transformation hyperpriors and the full-conditional distributions for γ γ γ are given in Appendix B.1.
Section 2.2 is flexible enough to allow for a transformation prior that implies cross-dependence among the elements of h, but we do not consider this case in this article. The main reason for this choice is that transformations are used in place of the original data set when implementing the preferred model (Step 4 of Algorithm 1). That is, the transformed values are used as a proxy for (or in place of) the data in the preferred model. Consequently, we would like to specify h to "overfit" the data so that h can be thought of as a proxy for the data. This leads to the following result, which shows that the prior in (6) leads to posterior replicates that overfit the data.
Proposition 3: Suppose (5) and (6) holds. Then, lim
. ., I 1 , j = 1, . . ., I 2 , and k = 1, . . ., I 3 .
Proof: See Appendix A.
Thus, the posterior mean of h (on the original scale of the data) are exactly the observed data {Z i j } as the hyperparameters go to zero. This suggests that estimates from the posterior distribution stated in Proposition 2 overfit the data. We emphasize that while f (h|z) overfits the data, it is not necessarily true that f (y, θ θ θ |z) overfits the data.
Assessment of the goodness of fit can be done similar to residual analyses of transformed data in traditional regression analyses. We compute the residuals δ δ δ = δ i j : i = 1, . . ., I j , j = 1, 2, 3 ′ , 
and where α is prespecified. A default choice is α = 0.05. In practice, it is rather straightforward to approximate q L,i j and q U,i j . Let h [b] i j and Y [b] i j be the b-th posterior replicate of h i j and Y i j so that δ
i j is the b-th posterior replicate of δ i j . Then q L,i j and q U,i j can be approximated with the α/2 and 1 − α/2 percentiles of the set {δ [b] i j : b = 1, . . ., B}, respectively. If the value of zero lies within this interval (e.g., q L,i j < 0 < q U,i j ) for many values of i and j, then this suggests that the model for y provides a reasonable fit to this data set.
Proposition 3, shows that the posterior mean of the transformation models overfits the data, which we motivated as a way to avoid oversmoothing estimates of y and θ θ θ in Algorithm 1. However, the fact that the transformation model overfits is also important from the point-of-view of diagnostics. In particular, in the goodness-of-fit literature, overfitted values are often used as a proxy for the data. For example, in log-linear models the most parsimonious reduced model that is not significantly different (in terms of the deviance or chi-square statistic) from the saturated model (an overfitted model) is used for statistical inference (e.g., see Agresti, 2007 , for a standard reference). This is exciting because it provides a new way to conduct classical residual analysis in a Bayesian multi-response context. In particular, in Section 3 we give an example of plotting the (posterior median) residuals versus a useful covariate not included in the analysis to assess whether or not it should be included in a model.
Friedman (1991) introduced a simulation design, which has become a useful benchmark study (e.g., see Chipman et al., 2010, among others) . Let h(x 1,i j , . . ., x 10,i j ) = 10sin(πx 1,i j x 2,i j ) + 20(x 3,i − 0.5) 2 + 10x 4,i j + 5x 5,i ; i = 1, . . ., I, j = 1, 2, 3,
which includes two non-linear terms, two linear terms, and a non-linear interaction. We consider the following specifications of the distributional assumptions associated with the data: 1,i3 , . . . , x 10,i3 ) ) , for i = 1, . . . , I j . Methods are compared using the root mean squared error (RMSE),
is estimated using Monte-Carlo integration using 2,000 iterations with a burn-in of 1,000. For each Bayesian method, we let g −1 j (h) be the pointwise posterior mean of g −1 j (h). We fit the preferred model using covariates x 1,i j , x 3,i j , x 4,i j , . . ., x 10,i j , and hence, we consider the case were an important covariate is not observed (i.e., {x 2,i j }) and several unneeded covariates are included (i.e., {x 6,i j , . . ., x 10,i j } are not present in (10)). The omissions of {x 2,i j } when implementing our method is a slight departure from the original setup in Friedman (1991) . However, we feel that it is more realistic to assume that not all covariates are observed in practice, and will be a helpful choice for illustration. We specify x k,i j ∼ Uniform(0, 1), where Uniform(0, 1) is a shorthand for the uniform distribution over the interval [0, 1] and k = 1, . . . , 10. The preferred models are SME and BART (and an extension), whose implementation are described in Appendix B.2 and B.3, respectively. In the implementation of each preferred method, we allow each response type to have different regression coefficients. Figure 1 : A violin plot of the RMSE (y-axis) by method (x-axis) over 20 independent replicates of the data. The data are simulated as described in Section 3.1. Each method is implemented using Algorithm 1, except the method "Saturated."
In this section, we evaluate the predictive performance of our Bayesian model with unknown transformations in the multi-response setting. In particular, we set the preferred model equal to BART (Chipman et al., 2010) and a Bayesian version of the SME model (Cressie and Johannesson, 2008) using basis functions introduced by (Hughes and Haran, 2013) . The posterior mean of h i j (referred to as the saturated model) are included as a default poor estimator, since it is known to overfit the data (see Proposition 3).
The data are simulated according to (11), with I = 1000, I 1 = 350, I 2 = 350, and I 3 = 200.
We repeat this simulation study 20 times, and we provide violin plots of the RMSE over the 20 replicates by method in Figure 1 . In Figure 2 we also plot the true function versus the estimated function for a single replicate data set. Figures 1 and 2 suggest that the transformation-based SME (BART) performs well in terms of predictive performance. For the replicate in Figure 2 the transformation-based SME (and BART) model had 97% (94%) of the point-wise credible intervals of the elements of δ δ δ containing zero. The patterns observed in Figure 1 mimic the goodness-of-fit diagnostics, which is notable because the goodness-of-fit diagnostics are data driven (and hence can be used in practice) while Figure 1 is based on the unknown truth. These results suggests that the Bayesian transformations can be used to obtain predictions in the non-Gaussian setting using two standard models, and also has a useful built-in goodness-of-fit diagnostic.
Now, suppose we have observed the values of {x 2,i j }, and recall these covariates are not included in the analysis. In Figure 3 , we plot the posterior median of the residuals versus the covariate {x 2,i j } across the indexes i and j for a single replicate of the data set. The plot clearly indicates a sinusoidal or possibly quadratic pattern, which suggests that this behavior is not captured in our model for y. We know this to be true because {x 2,i j } is not included in our implementation, but the data was generated using {x 2,i j }. This is an illustration of how our approach provides a Bayesian Figure 3 : We simulate a single replicate of {Y i j } according to Section 3.1. Then a SME model is implemented using the specifications in Appendix B.2. This plot displays the posterior median of {δ i j } (see Section 2.4) versus x 2,i j , which is not included in our implementation of the SME model. A systematic pattern in this plot suggests that including x 2,i j would improve our analysis of y.
analog to a graphical technique from classical regression analysis (i.e., systematic patterns in residuals from a multiple regression versus a covariate suggest that the covariate should be included in the analysis).
In this simulation study we compare the predictive performance our Bayesian transformation approach to predictions from the preferred model itself. A straightforward way to do this is to restrict ourselves to the continuous data-only setting, in which both modeling paradigms can be implemented. The data are simulated according to (11), with I 1 = 800, I = 1000, and I 2 = I 3 = 0.
We repeat this simulation study 20 times, and we provide violin plots of the RMSE over the 20 replicates by method in Figure 4 . In this section, we include an additional predictor: soft BART (SBART; Linero and Yang, 2018 , see Appendix B.3 for more details). We again see that the Bayesian transformation versions of BART and SME outperform the saturated model, with the SME model clearly outperforming BART. Additionally, the Bayesian transformation version of BART and SME perform only slightly better than or the same as their non-transformed counterparts. Here we see that SBART performs worse than the saturated model in terms of RMSE. The
Bayesian transformation version of SBART does not perform noticeably different than SBART in terms of RMSE. Thus, in the continuous only setting, if the preferred model performs well (or poorly) one should expect the Bayesian transformation approach to perform well (or poorly). Recall that we can use the goodness-of-fit approach in Section 2.4 to assess when a method performs poorly in practice. For example, for a single replicate data set, we found that the percent of credible intervals of the elements of δ δ δ that contain zero (by method) are as follows: 99.8% (SME), 77.4%
(BART), and 58.1% (SBART). This produces the same rankings of the method in terms of RMSE. where the cumulative costs have exceeded $1 billion. In total, these 246 weather events have lead to public costs that exceed $1.6 trillion. As such, it is extremely important that we have the foresight to determine which locations are the most sensitive to high costs due to natural disasters.
Thus, in this article, we consider jointly analyzing costs due to natural disasters using Bayesian transformations, where the preferred model is the MSTM (Bradley et al., 2015) .
To accurately forecast this spatio-temporal statistical process one needs to address several challenges. In particular, dynamics (i.e., how spatial fields change over time), non-stationarity (i.e., how costs at nearby time-points may not be similar), and non-separability (i.e., the possibility of complex interactions between locations and times). The MSTM is a model that allows for such patterns, however this model is not immediately useful for public health costs, since the data clearly shows zero inflation (e.g., see Sellers and Raim, 2016, among others) . See Figure 5 , for a histogram of public costs due to natural disasters in Florida counties from 2000 to 2018, where a large mode is present at zero costs.
Let Z i1 denote non-zero public costs due to natural disasters, Z i2 = 1 when the public costs are non-zero, and Z i2 = 0 when the public health costs are zero. There is dependence between these two response types, since Z i1 observed only when Z i2 = 1. Model Z i2 |h i2 independently as a
Bernoulli random variable with probability of success
Consider the following model for Z i j :
where f (Z i1 |h i1 , h i2 , Z i2 = 1) is a normal distribution with mean Y i1 and variance v > 0. The current non-Gaussian alternatives to the MSTM can not be used to model in (12), since they are defined only for a single response type; namely, Poisson (Bradley et al., 2018) and multinomial (Bradley et al., 2019c) data, respectively. The Bayesian transformation version of the MSTM, where the MSTM is defined by r = 74 basis functions, is intercept-only, and the first order nearest neighbor adjacency matrix is used. These 75 basis functions are chosen to be covariates made available by NCEI. Implementation of the MSTM is discussed in Appendix B.4. We provide a list of these covariates in Appendix C.
We analyze public costs due to natural disasters in Florida using data from 2000 to 2018.
The year 2018 is held out for cross-validation purposes. In Figure 6 we plot the predicted values from 2000 to 2017, and the corresponding data (both on the log-scale), and we see that the model produces small in sample error. This is supported by the number of credible intervals (as defined in Section 3.4) that contain zero, which is 100%. A naive baseline predictor to compare to is the previous years (i.e., 2017) public costs due to natural disasters. The RMSE (between predicted 2018 costs and observed costs) of this baseline predictor is 625.49, while the RMSE of the values of the proposed model is much smaller at 11.8996. Compare both values to the variance of the 2018 Florida natural disaster costs given by 1.5925 × 10 5 . Thus, we see that the proposed method provides reasonable out-of-sample error, and may be useful for estimating future costs due to natural disasters.
In this paper, we have provided several methodological contributions to Bayesian statistics. First, we have developed a general all-purpose Bayesian model to analyze multiple responses (e.g., continuous, Binomial counts, and Poisson counts). Our approach allows one to directly incorporate their preferred Bayesian model without adjusting the implementation of the preferred model. Second, a composite sampling method to implement a Bayesian model with unknown transformations has been introduced. Third, we developed a general Bayesian analog to the classical comparison between a saturated model and a reduced model. This results in the use of classical residual analysis for assessing goodness-of-fit in Bayesian models for multi-response data.
In our simulations, an illustration was given of non-linear functional analysis of multiple response types using BART as the preferred model. Additionally, an illustration was given of a joint spatial analysis of multiple response types using the SME model as the preferred model. These results suggest that the prediction error of our approach is small (in terms of RMSE), and we can develop multi-response versions of several different preferred models seamlessly. Additionally, data driven goodness-of-fit diagnostics were able to lead to the same conclusion as the RMSE criterion that is unobserved in practice. Finally, an illustration was given of a joint multivariate spatio-temporal analysis of zero-inflated continuous responses using the MSTM as the preferred model. Specifically, we analyze records of natural disaster costs in the state of Florida. The ability to accurately forecast the cost of natural disasters may help public policy professionals prepare.
There are several considerations for future methodological development. In particular, our choice of the transformation models in Section 2.3 implies that the data are independent. This is true even when the preferred model implies dependence among the transformed data values. Consequently, we are interested in extensions to our model that simultaneously allow for dependence among the data values and dependence among the transformed data values.
where f (z, h) = f (z|h) f (h|γ γ γ) f (γ γ γ) dγ γ γ and we have used the assumption of conditional independence between z and (y, θ θ θ ) given h. Then dividing by f (z) = f (z|h) f (h|γ γ γ) f (γ γ γ) dh dγ γ γ yields,
which is the desired result.
Proof of Proposition 2: Versions of this proof can be found in Diaconis and Ylvisaker (1979) and Bradley et al. (2019a) . The two distributions in (5) associated with j = 2 and j = 3 are members of the natural exponential family (Lehmann and Casella, 1998) , which are of the form,
where c i2 = b i and c i3 = 1. Upon multiplying by (6) we have:
which proves the result for j = 2 and j = 3. For j = 1,
which completes the results.
Proposition 3: In Equation (8) we see that
which converges to Z i1 as α 1 and κ 1 approach zero. The expectation of a beta distribution is well known (Casella and Berger, 2002) , which from (8) gives us
which converges to Z i2 /b i as α 2 and κ 2 approach zero. Similarly, the expectation of a gamma distribution is well known (Casella and Berger, 2002) , which from (8) gives us E {g(h i3 )|Z i3 , γ γ γ} = E(w 3 |Z i1 , γ γ γ) = α 3 + Z i3 κ 3 + 1 , which converges to Z i3 as α 3 and κ 3 approach zero.
in (B.1.1).
Consider the following mixed effects model for the transformed data (e.g., see Cressie and Johannesson, 2008, among others):
Data Model : h i j |β β β , η η η, ξ i j (·) ind ∼ Normal x ′ i j β β β + S ′ i j η η η + ξ i j , σ 2 ;
Process Model 1 : η η η|σ 2 η ∼ Normal 0 0 0 r , σ 2 η I r ;
Process Model 2 : ξ i j |σ 2 ξ ind ∼ Normal 0, σ 2 ξ ;
Parameter Model 1 : σ 2 ∼ IG (α v , β v ) ;
Parameter Model 2 : β β β ∼ Normal 0 0 0 p , σ 2 β I p ;
Parameter Model 3 : σ 2 ξ ∼ IG α ξ , β ξ ;
Parameter Model 4 : σ 2 η ∼ IG (α η , β η ) ; i = 1, . . . I j , j = 1, 2, 3, (B.2.1) where x i j is a p-dimensional vector of known covariates, I r is a r × r identity matrix, 0 0 0 r is an r-dimensional vector of zeros, α v = α η = α ξ = 1, β v = β η = β ξ = 1, σ 2 β = 100, and ξ ξ ξ = (ξ 11 , . . ., ξ I 3 3 ) ′ . The hyperparameters are chosen so that the prior is relatively "flat" and we find that our results are robust to these specifications. In Algorithm 1, we set Y i j = x ′ i j β β β +S ′ i j η η η +ξ i j and γ γ γ = (β β β ′ , σ 2 , σ 2 ξ , σ 2 η ) ′ .
The r-dimensional real-valued vector S i j is defined to be the Moran's I basis function (Hughes and Haran, 2013) . The Moran's I basis functions (Griffith, 2000 (Griffith, , 2002 (Griffith, , 2004 are motivated as a way to remove confounding between β β β and η η η, and allow for dimension reduction. The basis functions are derived from the Moran's I operator used in spatial statistics (Moran, 1950) . Specifically, basis functions are specified to be in the orthogonal column space associated with the hat matrix X(X ′ X) −1 X ′ , where the N × p matrix X = x i j : i = 1, . . . , I, j = 1, 2, 3 . Define the Moran's I distributions for variance parameters are well-known and are as follows:
Step 4 
Section 4 presents a bivariate example with a continuous response, and a binary response. The covariates used in Section 4 are listed below. All quoted terms are defined by NCEI.
• Intercept
• Interaction between an indicator for ℓ (i.e., variable) and ACS 2017 five year estimates of county-level population.
• Interaction between an indicator for ℓ and ACS 2017 five year estimates of county-level median income.
• Interaction between an indicator for ℓ and an indicator for the occurrence of an "astronomical low tide."
• Interaction between an indicator for ℓ and an indicator for the occurrence of a "coastal flood."
• Interaction between an indicator for ℓ and an indicator for the occurrence of a "cold/wind chill."
• Interaction between an indicator for ℓ and an indicator for the occurrence of a "debris flow."