{'fraction_non_alphanumeric': 0.042618095006520924, 'fraction_numerical': 0.013019739605207897, 'mean_word_length': 4.570056636296479, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Conditional Transformation Models

Recent developments in statistical regression methodology shift away from pure mean regression towards distributional regression models. One important strand thereof is that of conditional transformation models (CTMs). CTMs infer the entire conditional distribution directly by applying a transformation function to the response conditionally on a set of covariates towards a simple log-concave reference distribution. Thereby, CTMs allow not only variance, kurtosis or skewness but the complete conditional distribution to depend on the explanatory variables. We propose a Bayesian notion of conditional transformation models (BCTMs) focusing on exactly observed continuous responses, but also incorporating extensions to randomly censored and discrete responses. Rather than relying on Bernstein polynomials that have been considered in likelihood-based CTMs, we implement a spline-based parametrization for monotonic effects that are supplemented with smoothness priors. Furthermore, we are able to benefit from the Bayesian paradigm via easily obtainable credible intervals and other quantities without relying on large sample approximations. A simulation study demonstrates the competitiveness of our approach against its likelihood-based counterpart but also Bayesian additive models of location, scale and shape and Bayesian quantile regression. Two applications illustrate the versatility of BCTMs in problems involving real world data, again including the comparison with various types of competitors.


Introduction
Regression is omnipresent in many statistical applications and an ongoing field in recent research on statistical methods. While, in principle, interest always lies in describing the conditional distribution P Y |X=x of a response Y given a set of explanatory variables X with observed realisations x, most traditional approaches target the conditional expectation E(Y |X = x) as the only characteristic of interest (e.g. generalized linear or additive models; Nelder and Wedderburn, 1972;Hastie and Tibshirani, 1990). One way to abolish this often unwarranted simplification are generalized additive models for location scale and shape (GAMLSS; Rigby and Stasinopoulos, 2005) allowing for flexible relationships between the covariates and all parameters of the response distribution via flexible additive predictors. In this framework, the researcher can select from a diverse set of parametric distributions for discrete, continuous, mixed and multivariate response distributions (Klein, Kneib, Lang, Sohn et al., 2015;Klein, Kneib, Klasen and Lang, 2015). However, deciding on a parametric response distribution can quickly become a burden as it imposes strong assumptions on the model if not done with great care. One approach that does not entail a fixed parametric form is quantile regression introduced by Koenker and Bassett (1978). Since a distribution is satisfyingly captured by a dense grid of quantiles, each of these quantiles is modelled linearly or additively through covariates (Horowitz and Lee, 2005). Bayesian versions were developed e.g. by Yu and Moyeed (2001); Waldmann et al. (2013).
In contrast to all approaches mentioned so far, transformation models aim to infer the conditional distribution function directly. In an attempt to draw a bigger picture, we recapitulate a brief history of transformation models, while slightly focusing on Bayesian implementations. For a tour de force that concentrates on the frequentist perspective, see e.g. Hothorn et al. (2014); Möst (2015); Hothorn et al. (2018). Every transformation model entails a monotonically increasing transformation function h that acts on the response and is designed to reframe an unknown distribution P(Y ≤ y) in terms of the transformation, i.e. P(h(Y ) ≤ h(y)). The advent of parametric transformation models goes back to the Box-Cox model (Box and Cox, 1964) which ignited an area of active research that is still lit to this day. One approach that avoids strong assumptions on the parametric form of the transformation function was introduced by Cheng et al. (1995). It inspired plenty of models that share the estimation of a linear transformation function h(y|x) = h Y (y) − x β where the baseline transformation h Y (y) is estimated semiparametrically in conjunction with a linear, covariate-dependent shift x β. Prominent representatives are the proportional odds or the proportional hazards model. One of the first Bayesian transformation models was proposed by Pericchi (1981). Mallick and Walker (2003) model the transformation function h semiparametrically using (Bayesian) Bernstein polynomials and Pólya trees for the estimation of accelerated failure time (AFT) models among others. Song and Lu (2012) use Bayesian P-splines (Lang and Brezger, 2004) for transformation models with additive shift effects and a Gaussian reference distribution. Besides continuous responses, James et al. (2021) allow for discrete ordered and mixed discrete/continuous outcomes in conjunction with linear covariate effects.
Although very powerful in a lot of applications, transformation models of this type are considerably hindered by the additivity assumption on the scale of the transformation function h(y|x) = h Y (y) + h x (x) where the explanatory variable can only contribute a shift of the baseline transformation h Y and can therefore influence the conditional location parameter only. One modern example that includes linear interactions of covariates and gained a lot of attention is distribution regression (Chernozhukov et al., 2013). Here, in the context of counterfactuals, the conditional transformation function h(y|x) = h Y (y) − x β(y) is supplemented with varying-coefficient type interactions x β(y) where the varying coefficients β(y) are estimated on basis of P(Y ≤ y|X = x) = E(1(Y ≤ y)|X = x). This connection allows to account for heteroskedasticity or other patterns that vary with the covariates. An even more flexible variant comes with conditional transformation models (CTMs) as introduced by Hothorn et al. (2014) which share the same goal and aim to obtain an estimator for the whole conditional distribution function.
In this article, we propose the class of Bayesian conditional transformation models (BCTMs). BCTMs can be understood as a Bayesian interpretation of the likelihood-based CTMs of Hothorn et al. (2018) via the most likely transformation (MLT) model. All three models have in common that they target the direct estimation of the distribution function of a response Y conditional on a set of covariates X = x by means of estimating the conditional transformation function h(y|x). Yet, Bayesian inference based on Markov chain Monte Carlo (MCMC) simulations additionally allows us to obtain exact inferences on all quantities of interest without relying on large sample approximations or bootstrap procedures. This can be of particular value in scenarios with smaller samples where the parameters themselves are of secondary interest compared to complex transformations thereof. As the MLT model, the BCTM can be applied to discrete and continuous responses in the presence of random censoring, but additionally includes smoothness penalties for high-dimensional effects induced by the prior supporting stable function estimates.
Bayesian principles in connection with the modularity of implementation make it straightforward to expand BTCMs towards more complex prior structures enabling effect selection or different shrinkage properties for example. The idea of using monotonic P-splines for parametrizing transformation functions has been explored before (see e.g., Song and Lu, 2012;Tang et al., 2018), but our approach is innovative in a Bayesian setting with higher-dimensional interactions involving y, where curvatures are often complex and where control over the penalization mechanism can contribute to a better understanding of the model.
To summarize, we • introduce BCTMs as a new model class, • apply a B-spline basis in conjunction with reparametrized basis coefficients to impose monotonicity on the conditional transformation function in the y direction (opposed to the MLT model which relies on simple Bernstein polynomials), • supplement the unreparameterized vector of basis coefficients with a partially improper multivariate Gaussian prior that enforces smoothness towards a straight line both for monotonic and for unrestricted nonlinear effects, • develop Bayesian posterior estimation based on Hamiltonian Monte Carlo (HMC; Neal et al., 2011;Betancourt, 2017) using the highly-efficient No-U-Turn Sampler (NUTS, Hoffman and Gelman, 2014) for the vector of basis coefficients, • implement Bayesian model selection, • evaluate the distribution recovery ability and validity of credible intervals for BCTMs in different simulations and compare them to its main competitors, and • demonstrate different aspects and the practical relevance of BCTMs in applications on cholesterol levels from the Framingham heart study and on leukemia survival times.
The rest of the paper is structured as follows: Sec. 2 introduces BCTMs as a model class consisting of several building blocks including prior assumptions and theoretical properties.
Sec. 3 describes posterior estimation including Bayesian model selection. Sections 5 and 4 contain simulations and applications, respectively. Sec. 6 provides a brief review of our findings and proposes several directions for future research. The Supplement contains proofs of theoretical results as well as additional results for the simulations, an additional application on lung cancer survial times from the Veteran's Administration Lung Cancer Trial and further details.

Bayesian Conditional Transformation Models
In a CTM, the cumulative distribution function (CDF) of a response Y ∈ S ⊂ R conditional on a set of covariates X is specified via where the covariate-dependent function h(y|x) : S → R is assumed to be monotonically increasing in y to transform the response such that it follows a pre-specified reference distribution with continuous distribution function F Z : R → [0, 1]. The reference distribution is independent of x and does not contain any unknown parameters to be estimated. In this way, a CTM is characterized by the choice of the reference distribution F Z and a suitable parameterisation of the transformation function h(y|x) such that an estimate of the latter yields an estimate of a possibly complex conditional cumulative distribution function (cCDF)F Y |X . We will discuss both ingredients in more detail below. Naturally, distinctive characteristics of the transformation function such as monotonicity and smoothness are mirrored in F Y |X=x , which is why h(y|x) has to be modelled with great care.
Following Hothorn et al. (2014), we assume an additive decomposition on the scale of the transformation function into J partial transformation functions, i.e.
where h j (y|x), in the broadest sense, can be understood as response-covariate interactions that are monotone only in direction of y. To ensure identifiability, the partial transformation functions involving nonlinear terms are centered around zero, resulting in the additive decomposition h(y|x) = β 0 + J j=1 h j (y|x) with overall intercept β 0 of the conditional transformation function, which we will notationally suppress for most of what follows.
In light of (2), it is important to stress that additivity of the transformation function is assumed on the transformed scale, i.e. there is no explicit differentiation between signal and noise as in Gaussian regression models with separable error term. Hence, CTMs come with the benefit of a straightforward entry point to modelling all moments of the response distribution implicitly as functions of x. In the realm of CTMs, the flexibility of h j (y|x) constitutes the scope of the impact a covariate is admitted to have on the whole cCDF.
We assume that each of the J partial transformation functions h j (·|x), can be approximated by a linear combination of basis functions such that h j (y|x) = c j (y, x) γ j , where γ j is a vector of basis coefficients. Later we assume monotonicity of each partial transformation function in y, i.e. h j (y|x) = ∂h j (y|x) ∂y ≥ 0, which is sufficient but not necessary for an overall monotonic transformation function h(y|x). The complete transformation function and its derivative with respect to y are now given by with bases c(y, x) = (c 1 (y, x) , . . . , c J (y, x) ) and c (y, x) = (c 1 (y, x) , . . . , c J (y, x) ) for the transformation function and its derivative, respectively, and the stacked vector of all basis coefficients γ = (γ 1 , . . . , γ J ) .
In the following subsection, we introduce a generic and flexible joint basis for c j (y, x) that does not entail strict assumptions about the relationship between the moments of the response distribution and the respective covariates. Prior distributions, specific bases for the covariate effects, the choice of the references distribution F Z and a formal definition of our BCTM are covered in Sec. 2.2 to 2.6.

Generic conditional transformation functions
Let a j (y) and b j (x) denote vectors containing basis function evaluations B j1d 1 (y), d 1 = 1, . . . D 1 and B j2d 2 (x), d 2 = 1, . . . , D 2 for the response and the covariates, respectively, such that a j (y) = (B j11 (y), . . . , B j1D 1 (y)), and b j (x) = (B j21 (y), . . . , B j2D 2 (y)). Denoting by ⊗ the usual Kronecker product, we then obtain the most general form of partial transformation function in a BCTM as c j (y, x) = (a j (y) ⊗ b j (x) ) , leading to Essentially, the Kronecker product establishes a parametric interaction by forming pairwise products of the basis functions B j1d 1 (y) and B j2d 2 (x). The derivative with respect to y is therefore also a tensor product involving the differentiated basis functions B j1d 1 (y) = ∂B j1d 1 (y)/∂y. Specific restrictions on the two components of the tensor product lead to interesting special cases of the partial transformation function: • Setting a j (y) ≡ 1 leads to simple shift effects that only depend on the covariates.
• Setting b j (x) ≡ 1 yields an effect of only y that induces changes of the distributional shape (up to other effects).
• Linear effects a j (y) = (1, y) induce varying coefficient type effects where covariate effects linearly interact with the responses and the response takes the role of the interaction variable while the covariates are the effect modifiers.
Restricting our generic model to leads to a location-shift transformation model that comprises various earlier transformation models as special cases. In this case, only the location of the transformed response depends on the covariates via h x (x) and higher moments are captured unconditionally by the monotonic transformation h Y (y). This model type is parametrized by restricting the joint basis to c(y, x) = ((a(y) ⊗ 1) , (1 ⊗ b(x) ) ) = (a(y) , b(x) ) resulting in the shift transformation model c(y, x) γ = a(y) γ 1 + b(x) γ 2 .
We are relying on B-splines for the response dimension while various alternatives are available for the covariate dimension (see Sec. 2.3 for details). The choice of B-splines for representing a j (y) is mainly determined by the availability of suitable reparameterisations of the corresponding basis coefficients that ensure monotonicity along y in the tensor product for the partial response transformations and well-studied smoothness properties. More precisely, we follow Pya and Wood (2015) and is an identity matrix of size D 2 , Σ D 1 is a lower triangular matrix of size D 1 with Σ D 1 ,kl = 0 if k < l and Σ D 1 ,kl = 1 if k ≥ l, and the vectorβ j is Starting with a vector β j = (β j11 , . . . , β j1D 2 , β j21 , . . . β j2D 2 , . . . , β D 1 1 , . . . , β jD 1 D 2 ) ∈ R D of unconstrained parameters, these choices ensure that the vector of basis coefficients γ j is strictly increasing along the response dimension y which, in turn, implies a tensor product effect that is monotonically increasing along y. The complete model vectors of basis coefficients are then given by β = (β 1 , . . . , β J ) andβ = (β 1 , . . . ,β J ) , while the overall model matrix Σ is block diagonal with Σ 1 , . . . , Σ J as diagonal elements. We formalize the monotonicity of h(y|x) along the y dimension in the following theorem.
Theorem 2.1 (Monotonically increasing transformation function along y). Let h(·|x) : S → R be the transformation function (2) with basis representation (3) and partial transformation functions h j as in (4). Let furthermore γ j = Σ jβ j withβ j as in (6) and Σ j = Σ D 1 ⊗ I D 2 as defined above. Then, h(·|x) is monotonically increasing, that is for all y 1 , y 2 ∈ S with y 1 < y 2 we have h(y 1 |x) ≤ h(y 2 |x).
A proof of Theorem 2.1 can be found in the Supp. Part A.
In contrast to the MLT model introduced by Hothorn et al. (2018) that uses Bernstein polynomials as a basis for nonlinear effects, the BCTM is supplemented with a smoothnessinducing penalty through its prior and is therefore in principle less restrained regarding the number of model terms J and functional complexity in direction of the covariates.
Of course, other basis function representations than B-splines are immediately conceivable for a j (y). The main requirements for a suitable specification include the ability to incorporate monotonicity constraints, the analytical availability of the basis functions and their derivatives, and the numerically stable evaluation of these. While B-splines fulfill these requirements, investigating other choices and their properties is a promising avenue for future research.

Prior specifications
Overfitting of unregularized splines can be avoided in our Bayesian framework by enforcing smoothness and regularization through shrinkage priors. For the special case of B-splines, Bayesian P-splines assign multivariate Gaussian priors to the regression coefficient vectors.
We follow Kneib et al. (2019) and adopt this principle to tensor product terms such that the prior for the coefficient vector β j associated with one of the partial transformation functions h j in (4) is multivariate Gaussian with expectation zero and precision matrix where K j1 and K j2 are potentially rank deficient prior precision matrices of dimensions (D 1 × D 1 ) and (D 2 × D 2 ), respectively, controlling the type of smoothness required along the response and the covariate dimension, respectively. For the response dimension, we set Pya and Wood, 2015;Pya, 2010). The prior precision matrix K 2j shrinks in the direction of the respective covariate and the specific choice depends on the considered covariate effect of interect (see Sec. 2.3 for some examples). For monotonic nonlinear effects, the resulting penalty is quadratic in the (non-exponentiated) parameters β j . This corresponds to log differences in γ jd 1 d 2 for d 1 > 2, such that a first order random walk prior penalizes the squared differences between adjacent β jd 1 d 2 , resulting in shrinkage towards a straight line, similar to second order random walk penalties for univariate P-splines (Pya and Wood, 2015). The complete prior precision matrix K ≡ K(τ 2 , ω) is given as the block diagonal matrix with matrices K j as diagonal elements. We formalize the prior for γ j in the following proposition.
improper multivariate Gaussian prior with generalized inverse K − j of the prior precision matrix in (7). Assume furthermore for notational simplicity that D 2 = 1 such that γ jd1d2 ≡ γ jd1 . Then, the prior for γ j is given by Proof. The proof follows directly by applying the multivariate change of variable theorem twice to the transformation g 2 • g 1 : The amount of smoothness induced by the precision matrix (7) is controlled by the over-all smoothing variance τ 2 j > 0 and the weight parameter ω j ∈ [0, 1]. Following Kneib et al.
(2019), we assume a discrete prior for the latter which has the advantage that generalized determinants of K j can be pre-computed which considerably facilitates the numerically efficient implementation while still enabling anisotropic amounts of smoothness along the response and the covariate dimension. A uniform prior on a moderate number of equispaced values is used as a default for ω j . For the smoothing variance τ 2 j , we consider two alternatives: Standard inverse gamma (IG) priors τ 2 j ∼ IG(a j , b j ) where the hyperparameters are chosen among popular combinations such as a j = 1, b j = 0.001 to mimic a weakly informative setting, and scale-dependent (SD) hyperpriors as suggested in Klein and Kneib (2016). The latter results in a Weibull prior for τ 2 j with shape parameter 0.5 and scale parameter θ determined from a scaling criterion on expected effect sizes. We transfer this concept to partially monotonic tensor product effects where, to achieve numerical stability, it is important to control the variation of those parameters exponentiated in (6). More precisely, we consider the scaling criterion P (max d 1 =2,...,D 1 ,d 2 =1,...,D 2 |β jd 1 ,d 2 | ≤ c) = 1 − α with user-specified values for c and α. To determine the marginal prior distribution of τ 2 j required to evaluate the scaling criterion, we follow a simulation-based approach to marginalize out any additional hyperparameters. From the support of the exponential function, c = 3 and α = 0.01 are useful standards also used later in our empirical studies.
In a last step, we collect all model parameters in the vector with joint prior π ϑ (ϑ) which can be factorized into products of the individual priors. The coefficient β 0 denotes the intercept of the model, while β, τ 2 , ω are used to denote all basis coefficients, smoothing variances and anisotropy weights, respectively.

Bases for the covariate effects
We highlight special cases of bases b(x) relevant for our applications: • Linear effects. The basis for linear effects of covariates x 1 , . . . , x p collected in x is b j (x) = (x 1 , . . . , x p ) and we use a non-informative prior with K 2j = 0. This also applies to the overall intercept β 0 when centering the partial transformation functions.
• Random effects. Random effects (or frailties) are based on a grouping indicator g ∈ {1, . . . , G}. The resulting G-dimensional basis vector b j (g) has entry one if y belongs to group g and zero otherwise and we set K 2j = I G for i.i.d. random effects.
• Discrete spatial effects. Similar to random effects, a spatial effect of a discrete spatial variable s ∈ {1, . . . , S} is constructed as an indicator with entries in the S-dimensional basis vector b j (s) set to one if y belongs to region s and zero otherwise. We induce spatial smoothing in form of a Gaussian Markov random field (GMRF Rue and Held, 2005) prior. The precision matrix K 2j reflects the spatial orientation of the data, i.e. we define two regions as neighbours if they share a common border.
For an overview on possible bases for discrete (count) responses we refer the reader to Carlan and Kneib (2022).

Choice of the reference distribution
As already stated, it is the task of the conditional transformation function h(y|x) to transform the response values conditionally on the explanatory variables x such that they follow the reference distribution F Z . In that light, F Z plays a similar role as the inverse known link function in prominent model classes such as generalized linear models, but is less restrictive in the sense that the resulting conditional distribution does not have to be of known type. From a modelling perspective, it guarantees that the resulting estimated con-ditional density function integrates to one without requiring complex constraints. Note

Transformation densities
In this section, we introduce the conditional transformation densities f Y (y|β) given the vector of basis coefficients β (before the reparameterization). To emphasize that γ is a partially nonlinear reparameterization of β, we write γ(β).

Continuous responses
The density and log-density can easily be derived from equation (1) together with the parametrization of h and h in equation (4) such that where f Z denotes the density of the chosen reference distribution. In theory, for any absolute continuous response distribution F Y and reference distribution F Z with log-concave density f Z , there exists a unique, monotonically increasing transformation function h, such that F Y |X=x (·) = F Z (h(·|x)) (see Cor. 1 of Hothorn et al., 2018).
It is important to note that both for univariate and bivariate effects involving y, the part of the effect that belongs to the null space of K j consists of all location shifts and linear effects in y. In the context of BCTMs with the popular choice F Z = Φ, this means that the penalty shrinks towards the Gaussian location-scale family. In other words, the null space of the rank-deficient precision matrix consists of all Gaussian conditional distribution functions. This observation can be put to use when considered from the perspective of SD priors for the variances as described in Klein and Kneib (2016).
Discrete ordinal responses In case of discrete ordinal responses with a finite sample space where Y ∈ {y 1 , . . . , y K }, the corresponding conditional density function is given by For countably infinite sample spaces (as e.g. for count data) with Y ∈ {y 1 , y 2 , y 3 , . . .}, the density is given by Censored responses The Bayesian conditional transformation model incorporates all forms of random censoring. In the presence of censored observations, only the likelihood has to be adapted while the transformation function remains the same. The likelihood contributions for right-, left, and interval-censored continuous or discrete observations respectively are then given by 1 − F Z (c(y, x) γ(β)) for y ∈ (y, ∞) "right censored" It is also possible to adapt densities for truncated observations (Hothorn et al., 2018).

Formal definition of BCTMs
We end this section with a formal definition of BCTMs.
Definition 2.3 (BCTM). The quadruple ϑ, F Z , c, π ϑ (·) of unknown model parameters ϑ, a choice for the basis c, the reference distribution F Z and joint prior π ϑ is called Bayesian conditional transformation model (BCTM).

Posterior and estimation via MCMC
Assuming conditional independence the joint posterior is given by To obtain samples from (12) we use an MCMC sampler consisting of three alternating steps: Step 1. Sample from p(β|τ 2 , ω, y) using the NUTS.
Step 2. For j = 1, . . . , J, sample from p(τ 2 j |β j , y) using a Gibbs sampler in case of an IG prior or iteratively weighted least squares (IWLS) proposals in case of SD priors.
Step 3. For j = 1, . . . , J, sample ω j with a Gibbs step from its discrete full conditional.
The resulting MCMC samples can then be used to estimate the conditional distribution Updating the basis coefficients at Step 1. Basis coefficients are updated jointly by sampling from the log full conditional where the first term arises from one of the likelihoods described in Sec. 2.5 and the second term arises from the Gaussian prior). High dimensionality and strong dependencies among coefficients (stemming partly from the monotonicity constraints) aggravate sampling from the posterior distribution. This is further exacerbated by the mixed linear-nonlinear dependence of the transformation function onβ, rendering e.g. random-walk Metropolis algorithms slow and inefficient. One possible remedy lies in including gradient information as done by HMC. This, however, comes with the drawback that two additional tuning parameters (step size and number of leapfrog steps L) have to be set manually. To avoid this tricky task, we implement NUTS with dual averaging (Nesterov, 2009) that uses Hamiltonian principles for efficient exploration of the target distribution of β in an adaptive fashion The adaptive nature of NUTS enables a streamlined estimation process, effectively abolishing the need for costly preliminary tuning runs at the expense of some additional computation time per iteration which is owed mainly to the more sophisticated proposals.
The required gradient of the unnormalized log-posterior of the basis coefficients vector β for continuous responses is given by where C is a diagonal matrix with entries C dd = 1 ifβ d = β d , C dd = exp(β d ) otherwise, and similar expressions can straightforwardly be derived for discrete or censored responses.
Potentially flat parts of a fitted transformation function based on the reparameterization in Sec. 2.1 demand the parametersβ = exp(β) to be close to zero and thus the corresponding β to approach minus infinity. For NUTS, this does not result in overflow errors, but can lead to divergent transitions in the sampling path and NUTS trees with large tree depth as the different curvatures demand very different step sizes. Using a non-centered parametrization (Papaspiliopoulos et al., 2007) as a remedy is not feasible in a straightforward manner, because of the nonlinear transformation in the coefficient vectors. Instead, we found it helpful to increase the goal acceptance rate, forcing the sampler to take smaller steps, which is a small price to pay for the non-occurrence of divergencies. If the problem persists it is possible to drop unidentified (i.e. reparameterized coefficients that should be close to zero) in each iteration judging by the eigenvalues of the matrix square root of the Hessian of the posterior at (12) in an efficient way (Pya and Wood, 2015).
Furthermore, we resort to augmented precision matrices, e.g.
] + 10 −6 I to ensure positive definiteness and therefore a soft threshold for coefficient variances (Andrinopoulou et al., 2018). The NUTS warm-up phase can often be supported by standardizing each covariate or by rescaling them to (0,1). Both measures can facilitate mass matrix adaption. Regarding sampling efficiency, we found that using SD priors for the smoothing variances can decrease run times and improve the effective sample size.
Updating the smoothing variances at Step 2. When using an IG prior for the smoothing variances, they can be updated directly with a Gibbs step from the full conditional For the SD prior, updates can be implemented via IWLS proposals of log-variances following Klein and Kneib (2016).
Updating the weights at Step 3. The updates of the weights are straightforward using Gibbs sampling due to their discrete prior structure (Kneib et al., 2019).
Computational details While BCTMs are pretty robust regarding the choice of hyperparameters, varying them can improve computational speed and stabilize estimates that involve a monotonicity constraint. All results shown in Secs. 4, 5 were obtained with 4, 000 MCMC iterations with a NUTS warm-up phase of 2, 000 and a burn-in of 2, 000. Computations were carried out in R version 4.1.0 (R Core Team, 2020). To improve computing time, parts of the sampler were programmed using Rcpp (Eddelbuettel and Balamuta, 2017). The MASS matrix adaption scheme was adopted from adnuts (Monnahan and Kristensen, 2018).

Estimation of the cCDF
The resulting S MCMC samples can be used to estimate the cCDF F Y |X=x (y) asF Y |X=x,y (y) =

Model choice and variable selection
For model selection, we use the Watanabe-Akaike information criterion (WAIC Watanabe, 2010). It can be seen as approximation to computationally expensive cross validation (CV) and is conveniently computed from s = 1, . . . , S posterior samples. We validated WAIC against CV in some of our applications and found good agreements that support using information criteria as the basis for model choice and variable selection.
The WAIC overcomes certain limitations of the DIC (DIC; Spiegelhalter et al., 2002) such as its dependence on the posterior mean as a specific point estimate or the potential of observing negative effective parameter counts. It is given by WAIC = (−2l WAIC + 2p WAIC ), Var(log(f Y (y i |β))). In the regression literature, information criteria like the DIC and the WAIC are primarily used to discriminate between different types of response distributions and predictor specifications (e.g. Klein, Kneib, Lang, Sohn et al., 2015). In the holistic approach of BCTMs, the transformation function determines both the response distribution and the "predictor" which is why it is sufficient to use information criteria to compare different (partial) transformation function specifications that differ in flexibility and interaction structure. We also considered the deviance information criterion (DIC) as introduced by Spiegelhalter et al. (2002) which yielded similar results and is therefore omitted in the following.

Simulations
We conducted simulations to evaluate the empirical performance of BCTMs to recover the true data generating process compared to several competing methods from the literature (Sec. 4.1) and to provide valid uncertainty estimates by means of coverage rates (Sec. 4.2).
Simulation design For datasets of size n = 200, we generate two covariates as i.i.d.
realizations via x 1 ∼ U [0, 1] as well as x 2 ∼ U [−2, 2]. The response y is assumed to follow a heteroscedastic varying coefficient model (VCM) such that an appropriate CTM has to emulate a Gaussian location-scale model under the premises that the mean depends on the nonlinear varying coefficient (x 1 + 0.5) −1 for x 2 and that the variance is a nonlinear function of x 1 . To analyse the stability in the presence of noise variables, we consider six scenarios, where p = 0, . . . , 5 i.i.d. realizations from the standard uniform U [0, 1] with zero influence on the response are added. The complete vector of covariates is denoted by x p = (x 1 , x 2 , x 3 , . . . , x p+2 ) .
Benchmark methods For each of the resulting six scenarios, we fit • Lin. BCTM ϑ, Φ, ((1, y) ⊗ (1, x p )) , π ϑ (ϑ) : a restricted BCTM consisting of simple linear interactions • Lin. MLT: a linear MLT of the same type • Full BCTM ϑ, Φ, (a(y) ⊗(b(x 1 ) , . . . , b(x p+2 ) )) , π ϑ (ϑ) : a nonlinear BCTM consisting of nonlinear interactions with basis dimension of 10 in a and b • Full MLT: a nonlinear MLT of the same type with Bernstein polynomials of order 10, i.e. with joint basis (a Bs (y) ⊗ (b Bs (x 1 ) , . . . , b Bs,10 (x p+2 ) )) • Oracle BAMLSS: a Gaussian location-scale BAMLSS based on model (13), i.e. η µ = β 0 + x 2 · f (x 1 ) + p k=0 f (x 2+k ) and η σ 2 = β 0 + f (x 1 ) and • BAMLSS QR: a Bayesian semiparametric quantile regression specification with nonlinear effects of all explanatory variables Further details on the specifications are given in Supp. Tab. C.7. It is important to stress that Lin. BCTM/MLT and Oracle BAMLSS have in common that they are restricted by design to the true (Gaussian) distribution. In addition, the Oracle BAMLSS is the only model that is supplemented with the true predictor for the variance in all scenarios.
Yet, despite being linear in the covariates on the scale of the transformation function, the Lin. BCTM/MLT are nonlinear on the scale of the response. Since this information is in general not available, we also include the Full BCTM/MLT.
Performance measures As a first measure of performance, we computed the mean absolute deviation (MAD) of the estimates of F Y |X=x (y) from the true probabilities over a grid of y, x 1 and x 2 MAD(x 1 , x 2 ) = 1 n n i=1 |F Y |X=x (y i ) −F Y |X=x (y i )| based on 100 replications. Fig. 1 summarizes the empirical distributions of the minimum, median and maximum MAD for all models that provide estimates for the complete cCDF, i.e. all but the BAMLSS QR.
As a second performance measure, we computed conditional quantiles of the fitted response distribution corresponding to a sequence of probabilities α via numerical inversion. Fig. 2 shows the deviations of these from their true counterparts together with similar results obtained via QR BAMLSS which was used as a benchmark. Third, as a measure of accuracy that concentrates on tail features of the distribution, Fig. 3 shows the quantile score function QS α (F −1 (α), y) = 2(I(y < F −1 (α)) − α)(α − y) at α = 0.05 and α = 0.95, where I(A) = 1 if A is true, and zero otherwise (Gneiting, 2011). Last, to measure the overall forecast accuracy, we plot the decomposition of the continuous ranked probability score (CRPS;   according to the CRPS. In summary, BCTMs enable proper and reliable modelling of the complete cCDF and its quantiles, while avoiding restrictive assumptions on the shape of the distribution. BCTMs therefore provide a strong competitor in situations where neither the exact predictor specification nor the type of the response distribution can be derived from a priori considerations.

Coverage rates
To compare BCTM and MLT from a different perspective, we consider empirical coverage rates of pointwise 95% credible/confidence intervals in a simulation setting that concentrates on the estimation of nonlinear covariate effects.

Applications
We illustrate the versatility of BCTMs in three applications that differ with respect to the chosen reference distribution and transformation effect types. The first one highlights the applicability of the BCTM in the presence of highly skewed data (Sec. 5.1). The remaining two are BCTMs for (right-censored) survival data in form of a (non-)proportional hazards (NPH) model with random or spatial frailties (Sec. 5.2) and a partial (
• ϑ, Φ, (a(y) ⊗ b(age) , (year , sex )) , π ϑ (ϑ) : a full BCTM for age, where a(y) and b(age) contain an intercept, the tensor product is centered around zero and b(age) consists of a 10-dimensional B-splines basis leading to The default uses IG priors for both models but variants also employ the SD priors as competitors. Furthermore, we considered both models augmented by patient-specific i.i.d. distributional parameters (location, scale, degrees of freedom, and skewness) given by η k = β k,0 + x sex β k,sex + x age β k,age + x year β k,year , k = 1, . . . , K.
A variant thereof also contains the patient-specific i.i.d. random effects, see Supp. Tab. B.1 for full details on all model specifications.
Model selection All models are compared to each other using the DIC, WAIC and logscores in Supp. Tab. B.2. The log-scores are based on 10-fold CV. Overall, all criteria favour the tensor product spline BCTM with random effect over the GAMLSS specifically tailored to skewed responses. In general, the inclusion of random effects seems essential for obtaining realistic models while only smaller improvements result from the consideration of tensor products rather than VCMs. Replacing IG priors with SD priors does only yield a small performance improvement for the models without random effects. However, applying the SD prior results in noticeabe improvements in effectiveness and stability of the sampler,   Fig. 7 shows an estimated heat map that was obtained from the tensor product model assuming nonlinear covariate effects. While the general result is similar to the one in Fig. 6, we see a reversal of the trend towards right-skewness at age ≈ 55.

Leukemia survival
The second analysis considers acute myeloid leukemia survival of n = 1043 patients (Henderson et al., 2002) with 184 survival times being right-censored. In addition to the impact of the prognostic factors age, sex , white blood cell count (wbc) and the Townsend score (tpi ), indicating less affluent residential areas for higher values, we investigate spatial patterns in form of the indicator b(s) for 24 administrative regions in North West England.
In a first step, we fitted linear PH models ϑ, F MEV , (a(t) ⊗ 1) , x )) , π ϑ (ϑ) both without and with random effect for the administrative districts, and where F MEV denotes the CDF of the minimum extreme value distribution. Next, to account for spatial dependencies through a GMRF, we fit a ϑ, F MEV , (a(t) ⊗ 1) , (1 ⊗ (b(s) , x )) , π ϑ (ϑ) resulting in the spatial proportional hazards model As a last expansion, we fitted non-spatial and spatial NPH models for age, i.e. ϑ, F MEV , (a(t) ⊗ b(age)) , x )) , π ϑ (ϑ) and ϑ, F MEV , (a(t) ⊗ b(age)) , (1 ⊗ (b(s) , x )) , π ϑ (ϑ) , respectively. All model specifications are compared via the WAIC in Tab. 1 and respective posterior mean estimates of the log-negative harzard ratios are presented in Tab. 2. As a baseline check, Tab. 2 also includes esimates from the MLT for which however only the PH model with and without random effects for the districts and the NPH model (without random effects) can be estimated using the R packages tram (Hothorn, 2022) and tramME (Tamási and Hothorn, 2022). Details on the specifications can be found in Supp. Part B.2. Estimates of these models are similar to the ones of the corresponding BCTM.
Since overall the WAIC favours the spatial PH model (bctm spat), Fig. 8 shows the resulting estimated conditional survivor functions defined as for different Townsend scores (Panel A) accompanied by a depiction of the estimated spatial effect (Panel B). It confirms the findings of Tab. 2, indicating that affluency (lower tpi) is associated with higher survival at all times. The spatial effect of association to a district is associated with lower survival for higher values, and is therefore hinting on a lower mortality cluster in the northwest and on a high-risk "belt" running from northeast to southwest. Finally, Tab. B.5 of the Supplement shows the estimated posterior means of the log-negative hazard ratios (collected in γ 2 ), medians and credible intervals of the same model. Similar to the results in Zhou et al. (2020), we find that tpi , age and wbc are significant risk factors for surviving leukemia. bands for different quantiles of Townsend score (the lower, the more affluent) for females with age = 49 and wbc = 38.6(mean) for the first region. (B) Posterior mean spatial frailties (the larger, the higher mean mortality). Results are based on the spatial PH model ϑ, F MEV , (a(t) ⊗ 1) , (1 ⊗ (b(s) , x )) , π ϑ (ϑ) .

Summary and discussion
Our Bayesian treatment of CTMs based on MCMC is attractive for an assortment of reasons. Sampling-based inference provides posterior samples for coefficients of the conditional transformation function which can be transferred to samples of the cCDF, but also to samples of any quantity of interest that relies on the cCDF. It is straightforward, for example, to obtain point estimates and credible intervals without having to dive into asymptotics.
Furthermore, the Bayesian paradigm offers a natural way to impose smoothness penalties on the crucial nonlinear transformation functions, a feature likelihood-based competing methods (such as the MLT) are lacking in the software. In this way, the BCTM is able to resemble and even expand upon models ranging from simple to complex in settings with continuous, discrete and censored data without requiring strong assumptions.
In flat regions of the curve however, the reparameterization of the basis coefficients is such that the untransformed β may be weakly identified, resulting in potentially inefficient sample runs. This issue is explicitlty tackled in McKay Curtis and Ghosh (2011) who use a spike and slab prior directly on the basis coefficients that resulted in zero coefficients for flat regions, but only allows nonlinear monotonic covariate effects in a Gaussian setting. In our approach, we considered scale-dependent hyperpriors to counter mixing problems, but expanding such investigations to a wider scope with interactions is certainly an interesting field for future research.
Instead of specifying F Z a priori, it could also be interesting to include it as an additional free parameter in the estimation process. Among others, Linton et al. (2008); Politis (2013) describe the situation of a "model free" paradigm where the reference distribution is estimated without invoking any predetermined model (but by fully parameterizing the transformation function). This restriction is alleviated by the fact that in theory, arbitrarily complex distributions can be transformed to a basic reference distribution as long as the transformation function h is flexible enough. Abandoning the additivity assumption in h(y|x) in favor of e.g. tensor spline interactions however, can become computationally costly and numerically unstable due to the high dimensionality of the resulting basis, but can also be tackled by estimating the reference distribution in conjunction with a simpler structure in the transformation function. The idea of a free F Z was investigated in a Bayesian setting by Walker and Mallick (1999) and Mallick and Walker (2003) for example who use a Pólya tree prior for a series of (unconditional) semiparametric transformation models. Embedding it in the BCTM framework would result in a potentially very powerful addition to model flexibility.
Finally, our occupation with the BCTM for this article and beyond assured us that it represents an alluring modern competitor in the race to capture more and more aspects of the response distribution beyond the mean.

SUPPLEMENTARY MATERIAL
supplement.pdf This supplement contains the proof of Theorem 2.1 and further results for simulations and applications.
Code to reproduce the results from the applications is available on request.