# Most Likely Transformations

Torsten Hothorn
Universität Zürich
&Lisa Möst
Universität München
&Peter Bühlmann
ETH Zürich

Hothorn, Möst, and Bühlmann Most Likely Transformations Most Likely
Transformations We propose and study properties of maximum likelihood
estimators in the class of conditional transformation models. Based on a
suitable explicit parameterisation of the unconditional or conditional
transformation function, we establish a cascade of increasingly complex
transformation models that can be estimated, compared and analysed in the
maximum likelihood framework. Models for the unconditional or conditional
distribution function of any univariate response variable can be set-up and
estimated in the same theoretical and computational framework simply by
choosing an appropriate transformation function and parameterisation thereof.
The ability to evaluate the distribution function directly allows us to
estimate models based on the exact likelihood, especially in the presence of
random censoring or truncation. For discrete and continuous responses, we
establish the asymptotic normality of the proposed estimators. A reference
software implementation of maximum likelihood-based estimation for conditional
transformation models that allows the same flexibility as the theory developed
here was employed to illustrate the wide range of possible applications.
transformation model, distribution regression, conditional distribution
function, conditional quantile function, censoring, truncation transformation
model, distribution regression, conditional distribution function, conditional
quantile function, censoring, truncation Torsten Hothorn
Institut für Epidemiologie, Biostatistik und Prävention
Universität Zürich
Hirschengraben 84, CH-8001 Zürich, Switzerland
<EMAIL_ADDRESS>
Lisa Möst
Institut für Statistik
Ludwig-Maximilians-Universität München
Ludwigstraße 33, DE-80539 München, Germany
Peter Bühlmann
Seminar für Statistik
ETH Zürich
Rämistrasse 101, CH-8092 Zürich, Switzerland

††A shortened version of this technical report was accepted for publication in
the _Scandinavian Journal of Statistics_ 2017-06-19.

##  1 Introduction

In a broad sense, we can understand all statistical models as models of
distributions or certain characteristics thereof, especially the mean. All
distributions $\mathbb{P}_{Y}$ for at least ordered responses $Y$ can be
characterised by their distribution, quantile, density, odds, hazard or
cumulative hazard functions. In a fully parametric setting, all these
functions have been specified up to unknown parameters, and the ease of
interpretation can guide us in looking at the appropriate function. In the
semi- and non-parametric contexts, however, the question arises how we can
obtain an estimate of one of these functions without assuming much about their
shape. For the direct estimation of distribution functions, we deal with
monotonic functions in the unit interval, whereas for densities, we need to
make sure that the estimator integrates to one. The hazard function comes with
a positivity constraint, and monotonicity is required for the positive
cumulative hazard function. These computationally inconvenient restrictions
disappear completely only when the log-hazard function is estimated, and this
explains the plethora of research papers following this path. However, the
lack of any structure in the log-hazard function comes at a price. A too-
erratic behaviour of estimates of the log-hazard function has to be prevented
by some smoothness constraint; this makes classical likelihood inference
impossible. The novel characterisation and subsequent estimation of
distributions via their transformation function in a broad class of
transformation models that are developed in this paper can be interpreted as a
compromise between structure (monotonicity) and ease of parameterisation,
estimation and inference. This transformation approach to modelling and
estimation allows standard likelihood inference in a large class of models
that have so far commonly been dealt with by other inference procedures.

Since the introduction of transformation models based on non-linear
transformations of some response variable by Box and Cox (1964), this
attractive class of models has received much interest. In regression problems,
transformation models can be understood as models for the conditional
distribution function and are sometimes referred to as “distribution
regression”, in contrast to their “quantile regression” counterpart
(Chernozhukov _et al._ , 2013). Traditionally, the models were actively
studied and applied in the analysis of ordered categorical or censored
responses. Recently, transformation models for the direct estimation of
conditional distribution functions for arbitrary responses received interest
in the context of counterfactual distributions (Chernozhukov _et al._ , 2013),
probabilistic forecasting (Gneiting and Katzfuss, 2014), distribution and
quantile regression (Leorato and Peracchi, 2015; Rothe and Wied, 2013),
probabilistic index models (Thas _et al._ , 2012) and conditional
transformation models (Hothorn _et al._ , 2014). The core idea of any
transformation model is the application of a strictly monotonic transformation
function $h$ for the reformulation of an unknown distribution function
$\mathbb{P}(Y\leq y)$ as $\mathbb{P}(h(Y)\leq h(y))$, where the unknown
transformation function $h$ is estimated from the data. Transformation models
have received attention especially in situations where the likelihood contains
terms involving the conditional distribution function $\mathbb{P}(Y\leq
y\mid\text{\boldmath$X$}=\text{\boldmath$x$})=F_{Z}(h(y\mid\text{\boldmath$x$}))$
with inverse link function $F_{Z}$, most importantly for censored, truncated
and ordered categorical responses. For partially linear transformation models
with transformation function
$h(y\mid\text{\boldmath$x$})=h_{Y}(y)+h_{\text{\boldmath$x$}}(\text{\boldmath$x$})$,
much emphasis has been given to estimation procedures treating the baseline
transformation $h_{Y}$ (e.g. the log-cumulative baseline hazard function in
the Cox model) as a high-dimensional nuisance parameter. Prominent members of
these estimation procedures are the partial likelihood estimator (Cox, 1975)
and approaches influenced by the estimation equations introduced by Cheng _et
al._ (1995) and Chen _et al._ (2002). Once an estimate for the shift
$h_{\text{\boldmath$x$}}$ is obtained, the baseline transformation $h_{Y}$ is
then typically estimated by the non-parametric maximum likelihood estimator
(see, for example, Cheng _et al._ , 1997). An overview of the extensive
literature on the simultaneous non-parametric maximum likelihood estimation of
$h_{Y}$ and $h_{\text{\boldmath$x$}}$, i.e. estimation procedures not
requiring an explicit parameterisation of $h_{Y}$, for censored continuous
responses is given in Zeng and Lin (2007).

An explicit parameterisation of $h_{Y}$ is common in models of ordinal
responses (Tutz, 2012). For survival times, Kooperberg _et al._ (1995) and
Kooperberg and Clarkson (1997) introduced a cubic spline parameterisation of
the log-conditional hazard function with the possibility of response-varying
effects and estimated the corresponding models by maximum likelihood. Crowther
and Lambert (2014) followed up on this suggestion and used restricted cubic
splines. Many authors studied penalised likelihood approaches for spline
approximations of the baseline hazard function in a Cox model, e.g. Joly _et
al._ (1998); Cai and Betensky (2003); Sabanés Bové and Held (2013) or Ma _et
al._ (2014). Less frequently, the transformation function $h_{Y}$ was modelled
directly. Mallick and Walker (2003); Chang _et al._ (2005) and McLain and
Ghosh (2013) used Bernstein polynomials for $h_{Y}$, and Royston and Parmar
(2002) proposed a maximum likelihood approach using cubic splines for
modelling $h_{Y}$ and also time-varying effects. The connection between these
different transformation models is difficult to see because most authors
present their models in the relatively narrow contexts of survival or ordinal
data. The lack of a general understanding of transformation models made the
development of novel approaches in this model class burdensome. Hothorn _et
al._ (2014) decoupled the parameterisation of the conditional transformation
function $h(y\mid\text{\boldmath$x$})$ from the estimation procedure, and
showed that many interesting and novel models can be understood as
transformation models. The boosting-based optimisation of proper scoring
rules, however, was only developed for uncensored und right-censored (Möst and
Hothorn, 2015) observations in the absence of truncation and requires the
numerical approximation of the true target function. In a similar spirit,
Chernozhukov _et al._ (2013) applied the connection $\mathbb{P}(Y\leq
y\mid\text{\boldmath$X$}=\text{\boldmath$x$})=\mathbb{E}(\mathds{1}(Y\leq
y)\mid\text{\boldmath$X$}=\text{\boldmath$x$})$ for estimation in the
response-varying effects transformation model $\mathbb{P}(Y\leq
y\mid\text{\boldmath$X$}=\text{\boldmath$x$})=F_{Z}(h_{Y}(y)-\text{\boldmath$x$}^{\top}\text{\boldmath$\beta$}(y))$;
this approach can be traced back to Foresi and Peracchi (1995).

A drawback of all but the simplest transformation models is the lack of a
likelihood estimation procedure. Furthermore, although important connections
to other models have been known for some time (Doksum and Gasko, 1990), it is
often not easy to see how broad and powerful the class of transformation
models actually is. We address these issues and embed the estimation of
unconditional and conditional distribution functions of arbitrary univariate
random variables under all forms of random censoring and truncation into a
common theoretical and computational likelihood-based framework. In a
nutshell, we show in Section 2 that all distributions can be generated by a
strictly monotonic transformation of some absolute continuous random variable.
The likelihood function of the transformed variable can then be characterised
by this transformation function. The parameters of appropriate
parameterisations of the transformation function, and thus the parameters of
the conditional distribution function in which we are interested, can then be
estimated by maximum likelihood under simple linear constraints that allow
classical asymptotic likelihood inference, as will be shown in Section 3. Many
classical and contemporary models are introduced as special cases of this
framework. In particular, all transformation models sketched in this
introduction can be understood and estimated in this novel likelihood-based
framework. Extensions of classical and contemporary transformation models as
well as some novel models are derived from our unified theoretical framework
of transformation functions in Section 4, and their empirical performance is
illustrated and evaluated in Section 5.

##  2 The Likelihood of Transformations

Let $(\Omega,\mathfrak{A},\mathbb{P})$ denote a probability space and
$(\Xi,\mathfrak{C})$ a measureable space with at least ordered sample space
$\Xi$. We are interested in inference about the distribution $\mathbb{P}_{Y}$
of a random variable $Y$, i.e. the probability space
$(\Xi,\mathfrak{C},\mathbb{P}_{Y})$ defined by the $\mathfrak{A}-\mathfrak{C}$
measureable function $Y:\Omega\rightarrow\Xi$. For the sake of notational
simplicity, we present our results for the unconditional and ordered case
first; regression models and unordered responses are discussed in Section 4.2.
The distribution $\mathbb{P}_{Y}=f_{Y}\odot\mu$ is dominated by some measure
$\mu$ and characterised by its density function $f_{Y}$, distribution function
$F_{Y}(y)=\mathbb{P}_{Y}(\\{\xi\in\Xi\mid\xi\leq y\\})$, quantile function
$F^{-1}_{Y}(p)=\inf\\{y\in\Xi\mid F_{Y}(y)\geq p\\}$, odds function
$O_{Y}(y)=F_{Y}(y)/(1-F_{Y}(y))$, hazard function
$\lambda_{Y}(y)=f_{Y}(y)/(1-F_{Y}(y))$ or cumulative hazard function
$\Lambda_{Y}(y)=-\log(1-F_{Y}(y))$. For notational convenience, we assume
strict monotonicity of $F_{Y}$, i.e. $F_{Y}(y_{1})<F_{Y}(y_{2})\,\forall
y_{1}<y_{2}\in\Xi$. Our aim is to obtain an estimate $\hat{F}_{Y,N}$ of the
distribution function $F_{Y}$ from a random sample
$Y_{1},\dots,Y_{N}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\mathbb{P}_{Y}$.
In the following, we will show that one can always write this potentially
complex distribution function $F_{Y}$ as the composition of a much simpler and
a priori specified distribution function $F_{Z}$ and a strictly monotonic
transformation function $h$. The task of estimating $F_{Y}$ is then reduced to
obtaining an estimate $\hat{h}_{N}$. The latter exercise, as we will show in
this paper, is technically and conceptually attractive.

Let $(\mathbb{R},\mathfrak{B})$ denote the Euclidian space with Borel
$\sigma$-algebra and $Z:\Omega\rightarrow\mathbb{R}$ an
$\mathfrak{A}-\mathfrak{B}$ measureable function such that the distribution
$\mathbb{P}_{Z}=f_{Z}\odot\mu_{L}$ is absolutely continuous ($\mu_{L}$ denotes
the Lebesgue measure) in the probability space
$(\mathbb{R},\mathfrak{B},\mathbb{P}_{Z})$. Let $F_{Z}$ and $F^{-1}_{Z}$
denote the corresponding distribution and quantile functions. We furthermore
assume $0<f_{Z}(z)<\infty\,\forall z\in\mathbb{R}$, $F_{Z}(-\infty)=0$ and
$F_{Z}(\infty)=1$ for a log-concave density $f_{Z}$ as well as the existence
of the first two derivatives of the density $f_{Z}(z)$ with respect to $z$;
both derivatives shall be bounded. We do not allow any unknown parameters for
this distribution. Possible choices include the standard normal, standard
logistic (SL) and minimum extreme value (MEV) distribution with distribution
functions $F_{Z}(z)=\Phi(z)$,
$F_{Z}(z)=F_{\operatorname{SL}}(z)=(1+\exp(-z))^{-1}$ and
$F_{Z}(z)=F_{\operatorname{MEV}}(z)=1-\exp(-\exp(z))$, respectively. In the
first step, we will show that there always exists a unique and strictly
monotonic transformation $g$ such that the unknown and potentially complex
distribution $\mathbb{P}_{Y}$ that we are interested in can be generated from
the simple and known distribution $\mathbb{P}_{Z}$ via
$\mathbb{P}_{Y}=\mathbb{P}_{g\circ Z}$. More formally, let
$g:\mathbb{R}\rightarrow\Xi$ denote a $\mathfrak{B}-\mathfrak{C}$ measureable
function. The composition $g\circ Z$ is a random variable on
$(\Xi,\mathfrak{C},\mathbb{P}_{g\circ Z})$. We can now formulate the existence
and uniqueness of $g$ as a corollary to the probability integral transform.

######  Corollary 1.

For all random variables $Y$ and $Z$, there exists a unique strictly
monotonically increasing transformation $g$, such that
$\mathbb{P}_{Y}=\mathbb{P}_{g\circ Z}$.

###### Proof.

Let $g=F_{Y}^{-1}\circ F_{Z}$ and $Z\sim\mathbb{P}_{Z}$. Then
$U:=F_{Z}(Z)\sim\operatorname{U}[0,1]$ and $Y=F_{Y}^{-1}(U)\sim\mathbb{P}_{Y}$
by the probability integral transform. Let $h:\Xi\rightarrow\mathbb{R}$, such
that $F_{Y}(y)=F_{Z}(h(y))$. From

$\displaystyle F_{Y}(y)=(F_{Z}\circ F_{Z}^{-1}\circ F_{Y})(y)\iff
h=F_{Z}^{-1}\circ F_{Y}$

we get the uniqueness of $h$ and therefore $g$. The quantile function
$F_{Z}^{-1}$ and the distribution function $F_{Y}$ exist by assumption and are
both strictly monotonic and right-continuous. Therefore, $h$ is strictly
monotonic and right-continuous and so is $g$. ∎

######  Corollary 2.

For $\mu=\mu_{L}$, we have $g=h^{-1}$ and $h^{\prime}(y)=\frac{\partial
h(y)}{\partial y}=f_{Z}((F_{Z}^{-1}\circ F_{Y})(y))^{-1}f_{Y}(y)$.

This result for absolutely continuous random variables $Y$ can be found in
many textbooks (e.g. Lindsey, 1996), Corollary 1 also covers the discrete
case.

######  Corollary 3.

For the counting measure $\mu=\mu_{C}$, $h=F_{Z}^{-1}\circ F_{Y}$ is a right-
continuous step function because $F_{Y}$ is a right-continuous step function
with steps at $y\in\Xi$.

We now characterise the distribution $F_{Y}$ by the corresponding
transformation function $h$, set up the corresponding likelihood of such a
transformation function and estimate the transformation function based on this
likelihood. Let
$\mathcal{H}=\\{h:\Xi\rightarrow\mathbb{R}\mid\mathfrak{C}-\mathfrak{B}\text{
measureable},h(y_{1})<h(y_{2})\,\forall y_{1}<y_{2}\in\Xi\\}$ denote the space
of all strictly monotonic transformation functions. With the transformation
function $h$, we can evaluate $F_{Y}$ as $F_{Y}(y\mid h)=F_{Z}(h(y))\,\forall
y\in\Xi$. Therefore, we only need to study the transformation function $h$;
the inverse transformation $g=h^{-1}$ (used to define a “group family” or
“group model” by Lehmann, 1983; Bickel _et al._ , 1993) is not necessary in
what follows. The density for absolutely continuous variables $Y$
($\mu=\mu_{L}$) is now given by

$\displaystyle f_{Y}(y\mid h)=f_{Z}(h(y))h^{\prime}(y).$

For discrete responses $Y$ ($\mu=\mu_{C}$) with finite sample space
$\Xi=\\{y_{1},\dots,y_{K}\\}$, the density is

$\displaystyle f_{Y}(y_{k}\mid
h)=\left\\{\begin{array}[]{ll}F_{Z}(h(y_{k}))&k=1\\\
F_{Z}(h(y_{k}))-F_{Z}(h(y_{k-1}))&k=2,\dots,K-1\\\
1-F_{Z}(h(y_{k-1}))&k=K\end{array}\right.$

and for countably infinite sample spaces $\Xi=\\{y_{1},y_{2},y_{3},\dots\\}$,
we get the density

$\displaystyle f_{Y}(y_{k}\mid
h)=\left\\{\begin{array}[]{ll}F_{Z}(h(y_{k}))&k=1\\\
F_{Z}(h(y_{k}))-F_{Z}(h(y_{k-1}))&k>1.\end{array}\right.$

With the conventions $F_{Z}(h(y_{0})):=F_{Z}(h(-\infty)):=0$ and
$F_{Z}(h(y_{K})):=F_{Z}(h(\infty)):=1$, we use the more compact notation
$f_{Y}(y_{k}\mid h)=F_{Z}(h(y_{k}))-F_{Z}(h(y_{k-1}))$ in the sequel.

For a given transformation function $h$, the likelihood contribution of a
datum $C=(\underaccent{\bar}{\ry},\bar{y}]\in\mathfrak{C}$ is defined in terms
of the distribution function (Lindsey, 1996):

$\displaystyle\mathcal{L}(h\mid Y\in C):=\int_{C}f_{Y}(y\mid
h)d\mu(y)=F_{Z}(h(\bar{y}))-F_{Z}(h(\underaccent{\bar}{\ry})).$

This “exact” definition of the likelihood applies to most practical situations
of interest and, in particular, allows discrete and (conceptually) continuous
as well as censored or truncated observations $C$. For a discrete response
$y_{k}$, we have $\bar{y}=y_{k}$ and $\underaccent{\bar}{\ry}=y_{k-1}$, such
that $\mathcal{L}(h\mid Y=y_{k})=f_{Y}(y_{k}\mid
h)=F_{Z}(h(\bar{y}))-F_{Z}(h(\underaccent{\bar}{\ry}))$. For absolutely
continuous random variables $Y$, we almost always observe an imprecise datum
$(\underaccent{\bar}{\ry},\bar{y}]\subset\mathbb{R}$ and, for short intervals
$(\underaccent{\bar}{\ry},\bar{y}]$, approximate the exact likelihood
$\mathcal{L}(h\mid Y\in(\underaccent{\bar}{\ry},\bar{y}])$ by the term
$(\bar{y}-\underaccent{\bar}{\ry})f_{Y}(y\mid h)$ or simply $f_{Y}(y\mid h)$
with $y=(\underaccent{\bar}{\ry}+\bar{y})/2$ (Lindsey, 1999). This
approximation only works for relatively precise measurements, i.e. short
intervals. If longer intervals are observed, one speaks of “censoring” and
relies on the exact definition of the likelihood contribution instead of using
the above approximation (Klein and Moeschberger, 2003). In summary, the
likelihood contribution of a conceptually “exact continuous” or left-, right-
or interval-censored continuous or discrete observation
$(\underaccent{\bar}{\ry},\bar{y}]$ is given by

$\displaystyle\mathcal{L}(h\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])\left\\{\begin{array}[]{ll}\approx
f_{Z}(h(y))h^{\prime}(y)&y=(\underaccent{\bar}{\ry}+\bar{y})/2\in\Xi\quad\text{```exact
continuous'''}\\\
=1-F_{Z}(h(\underaccent{\bar}{\ry}))&y\in(\underaccent{\bar}{\ry},\infty)\cap\Xi\quad\text{`right-
censored'}\\\ =F_{Z}(h(\bar{y}))&y\in(-\infty,\bar{y}]\cap\Xi\quad\text{`left-
censored'}\\\
=F_{Z}(h(\bar{y}))-F_{Z}(h(\underaccent{\bar}{\ry}))&y\in(\underaccent{\bar}{\ry},\bar{y}]\cap\Xi\quad\text{`interval-
censored',}\end{array}\right.$

under the assumption of random censoring. The likelihood is more complex under
dependent censoring (Klein and Moeschberger, 2003), but we will not elaborate
on this issue. The likelihood contribution $\mathcal{L}(h\mid
Y\in(y_{k},y_{k-1}])$ of an ordered factor in category $y_{k}$ is equivalent
to the term $\mathcal{L}(h\mid Y\in(\underaccent{\bar}{\ry},\bar{y}])$
contributed by an interval-censored observation
$(\underaccent{\bar}{\ry},\bar{y}]$, when category $y_{k}$ is defined by the
interval $(\underaccent{\bar}{\ry},\bar{y}]$. Thus, the expression
$F_{Z}(h(\bar{y}))-F_{Z}(h(\underaccent{\bar}{\ry}))$ for the likelihood
contribution reflects the equivalence of interval-censoring and categorisation
at corresponding cut-off points.

For truncated observations in the interval $(y_{l},y_{r}]\subset\Xi$, the
above likelihood contribution is defined in terms of the distribution function
conditional on the truncation

$\displaystyle F_{Y}(y\mid Y\in(y_{l},y_{r}])=F_{Z}(h(y)\mid
Y\in(y_{l},y_{r}])=\frac{F_{Z}(h(y))}{F_{Z}(h(y_{r}))-F_{Z}(h(y_{l}))}\quad\forall
y\in(y_{l},y_{r}]$

and thus the likelihood contribution changes to (Klein and Moeschberger, 2003)

$\displaystyle\frac{\mathcal{L}(h\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])}{F_{Z}(h(y_{r}))-F_{Z}(h(y_{l}))}=\frac{\mathcal{L}(h\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])}{\mathcal{L}(h\mid
Y\in(y_{l},y_{r}])}\quad\text{when }y_{l}<\underaccent{\bar}{\ry}<\bar{y}\leq
y_{r}.$

It is important to note that the likelihood is always defined in terms of a
distribution function (Lindsey, 1999), and it therefore makes sense to
directly model the distribution function of interest. The ability to uniquely
characterise this distribution function by the transformation function $h$
gives rise to the following definition of an estimator $\hat{h}_{N}$.

######  Definition 1 (Most likely transformation).

Let $C_{1},\dots,C_{N}$ denote an independent sample of possibly randomly
censored or truncated observations from $\mathbb{P}_{Y}$. The estimator

$\displaystyle\hat{h}_{N}:=\operatorname*{{arg\,max}}_{\tilde{h}\in\mathcal{H}}\sum_{i=1}^{N}\log(\mathcal{L}(\tilde{h}\mid
Y\in C_{i}))$

is called the most likely transformation (MLT).

Log-concavity of $f_{Z}$ ensures concavity of the log-likelihood (except when
all observations are right-censored) and thus ensures the existence and
uniqueness of $\hat{h}_{N}$.

##### Example

For an absolutely continuous response $Y$ the likelihood and log-likelihood
for $h$ are approximated by the density and log-density evaluated at
$y=(\underaccent{\bar}{\ry}+\bar{y})/2$, respectively:

$\displaystyle\mathcal{L}(h\mid Y\in(\underaccent{\bar}{\ry},\bar{y}])$
$\displaystyle\approx$ $\displaystyle f_{Z}(h(y))h^{\prime}(y)$
$\displaystyle\log(\mathcal{L}(h\mid Y\in(\underaccent{\bar}{\ry},\bar{y}]))$
$\displaystyle\approx$ $\displaystyle\log(f_{Z}(h(y)))+\log(h^{\prime}(y)).$

Strict monotonicity of the transformation function $h$ is required; otherwise
the likelihood is not defined. The term $\log(h^{\prime}(y))$ is not a penalty
term, but the likelihood favours transformation functions with a large
positive derivative at the observations. If we assume
$Y\sim\operatorname{N}(\alpha,\sigma^{2})$ and for the choice
$Z\sim\operatorname{N}(0,1)$ with $F_{Z}=\Phi$ and $f_{Z}=\phi$, we can
restrict $h$ to linear functions $h(y)=(y-\alpha)\sigma^{-1}$. The likelihood
reduces to

$\displaystyle\mathcal{L}(h\mid Y\in(\underaccent{\bar}{\ry},\bar{y}])\approx
f_{Z}(h(y))h^{\prime}(y)=\phi((y-\alpha)\sigma^{-1})\sigma^{-1}=\phi_{\alpha,\sigma^{2}}(y)=f_{Y}(y\mid\alpha,\sigma^{2}).$

In this simple location-scale family, the most likely transformation is
characterised by the parameters of the normal distribution of $Y$. It is
important to note that for other choices of $F_{Z}$, the most likely
transformation is non-linear; however, the distribution function
$F_{Y}=F_{Z}(h(y))$ is invariant with respect to $F_{Z}$ because we can always
write $h$ as $F_{Z}^{-1}\circ F_{Y}$. In other words, with $F_{Z}\neq\Phi$, we
can still model normal responses $Y$; however, a non-linear transformation
function $h$ is required.

Many distributions are defined by a transformation function $h$, for example,
the Box-Cox power exponential family (Rigby and Stasinopoulos, 2004), the
sinh-arcsinh distributions (Jones and Pewsey, 2009), or the T-X family of
distributions (Alzaatreh _et al._ , 2013). The parameters of these
distributions can, for example, be estimated by the GAMLSS approach (Rigby and
Stasinopoulos, 2005). In what follows, we do not assume any specific form of
the transformation function but parameterise $h$ in terms of basis functions.
We now introduce such a parameterisation, a corresponding family of
distributions, a maximum likelihood estimator and a large class of models for
unconditional and conditional distributions.

##  3 Transformation Analysis

We parameterise the transformation function $h(y)$ as a linear function of its
basis-transformed argument $y$ using a basis function
$\text{\boldmath$a$}:\Xi\rightarrow\mathbb{R}^{P}$, such that
$h(y)=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$},\text{\boldmath$\vartheta$}\in\mathbb{R}^{P}$.
The choice of the basis function $a$ is problem specific and will be discussed
in Section 4. The likelihood $\mathcal{L}$ only requires evaluation of $h$,
and only an approximation thereof using the Lebesgue density of “exact
continuous” observations makes the evaluation of the first derivative of
$h(y)$ with respect to $y$ necessary. In this case, the derivative with
respect to $y$ is given by
$h^{\prime}(y)=\text{\boldmath$a$}^{\prime}(y)^{\top}\text{\boldmath$\vartheta$}$,
and we assume that $\text{\boldmath$a$}^{\prime}$ is available. In the
following, we will write
$h=\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$ and
$h^{\prime}={\text{\boldmath$a$}^{\prime}}^{\top}\text{\boldmath$\vartheta$}$
for the transformation function and its first derivative, omitting the
argument $y$, and we assume that both functions are bounded away from
$-\infty$ and $\infty$. For a specific choice of $F_{Z}$ and $a$, the
transformation family of distributions consists of all distributions
$\mathbb{P}_{Y}$ whose distribution function $F_{Y}$ is given as the
composition $F_{Z}\circ\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$;
this family can be formally defined as follows.

######  Definition 2 (Transformation family).

The distribution family

$\displaystyle\mathbb{P}_{Y,\Theta}=\\{F_{Z}\circ\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}\mid\text{\boldmath$\vartheta$}\in\Theta\\}$

with parameter space
$\Theta=\\{\text{\boldmath$\vartheta$}\in\mathbb{R}^{P}\mid\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}\in\mathcal{H}\\}$
is called transformation family of distributions
$\mathbb{P}_{Y,\text{\boldmath$\vartheta$}}$ with transformation functions
$\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}\in\mathcal{H}$,
$\mu$-densities $f_{Y}(y\mid\text{\boldmath$\vartheta$}),y\in\Xi$, and error
distribution function $F_{Z}$.

The classical definition of a transformation family relies on the idea of
invariant distributions, i.e. only the parameters of a distribution are
changed by a transformation function but the distribution itself is not
changed. The normal family characterised by affine transformations is the most
well-known example (e.g. Fraser, 1968; Lindsey, 1996). Here, we explicitly
allow and encourage transformation functions that change the shape of the
distribution. The transformation function
$\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$ is, at least in
principle, flexible enough to generate any distribution function
$F_{Y}=F_{Z}\circ\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$ from
the distribution function $F_{Z}$. We borrow the term “error distribution”
function for $F_{Z}$ from Fraser (1968), because $Z$ can be understood as an
error term in some of the models discussed in Section 4. The problem of
estimating the unknown transformation function $h$, and thus the unknown
distribution function $F_{Y}$, reduces to the problem of estimating the
parameter vector $\vartheta$ through maximisation of the likelihood function.
We assume that the basis function $a$ is such that the parameters $\vartheta$
are identifiable.

######  Definition 3 (Maximum likelihood estimator).

$\displaystyle\hat{\text{\boldmath$\vartheta$}}_{N}:=\operatorname*{{arg\,max}}_{\text{\boldmath$\vartheta$}\in\Theta}\sum_{i=1}^{N}\log(\mathcal{L}(\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}\mid
Y\in C_{i}))$

Based on the maximum likelihood estimator
$\hat{\text{\boldmath$\vartheta$}}_{N}$, we define plug-in estimators of the
most likely transformation function and the corresponding estimator of our
target distribution $F_{Y}$ as
$\hat{h}_{N}:=\text{\boldmath$a$}^{\top}\hat{\text{\boldmath$\vartheta$}}_{N}$
and $\hat{F}_{Y,N}:=F_{Z}\circ\hat{h}_{N}$. Because the problem of estimating
an unknown distribution function is now embedded in the maximum likelihood
framework, the asymptotic analysis benefits from standard results on the
asymptotic behaviour of maximum likelihood estimators. We begin with deriving
the score function and Fisher information. The score contribution of an “exact
continuous” observation $y=(\underaccent{\bar}{\ry}+\bar{y})/2$ from an
absolutely continuous distribution is approximated by the gradient of the log-
density

$\displaystyle\text{\boldmath$s$}(\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])\approx\frac{\partial\log(f_{Y}(y\mid\text{\boldmath$\vartheta$}))}{\partial\text{\boldmath$\vartheta$}}$
$\displaystyle=$
$\displaystyle\frac{\partial\log(f_{Z}(\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$})))+\log({\text{\boldmath$a$}^{\prime}(y)}^{\top}\text{\boldmath$\vartheta$})}{\partial\text{\boldmath$\vartheta$}}$
(4) $\displaystyle=$
$\displaystyle\text{\boldmath$a$}(y)\frac{f_{Z}^{\prime}(\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$})}{f_{Z}(\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$})}+\frac{\text{\boldmath$a$}^{\prime}(y)}{{\text{\boldmath$a$}^{\prime}(y)}^{\top}\text{\boldmath$\vartheta$}}.$

For an interval-censored or discrete observation $\underaccent{\bar}{\ry}$ and
$\bar{y}$ (the constant terms
$F_{Z}(\text{\boldmath$a$}(-\infty)^{\top}\text{\boldmath$\vartheta$})=F_{Z}(-\infty)=0$
and
$F_{Z}(\text{\boldmath$a$}(\infty)^{\top}\text{\boldmath$\vartheta$})=F_{Z}(\infty)=1$
vanish), the score contribution is

$\displaystyle\text{\boldmath$s$}(\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])$ $\displaystyle=$
$\displaystyle\frac{\partial\log(\mathcal{L}(\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}]))}{\partial\text{\boldmath$\vartheta$}}$
(5) $\displaystyle=$
$\displaystyle\frac{\partial\log(F_{Z}(\text{\boldmath$a$}(\bar{y})^{\top}\text{\boldmath$\vartheta$})-F_{Z}(\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}\text{\boldmath$\vartheta$}))}{\partial\text{\boldmath$\vartheta$}}$
$\displaystyle=$
$\displaystyle\frac{f_{Z}(\text{\boldmath$a$}(\bar{y})^{\top}\text{\boldmath$\vartheta$})\text{\boldmath$a$}(\bar{y})-f_{Z}(\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}\text{\boldmath$\vartheta$})\text{\boldmath$a$}(\underaccent{\bar}{\ry})}{F_{Z}(\text{\boldmath$a$}(\bar{y})^{\top}\text{\boldmath$\vartheta$})-F_{Z}(\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}\text{\boldmath$\vartheta$})}.$

For a truncated observation, the score function is
$\text{\boldmath$s$}(\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])-\text{\boldmath$s$}(\text{\boldmath$\vartheta$}\mid
Y\in(y_{l},y_{r}])$.

The contribution of an “exact continuous” observation $y$ from an absolutely
continuous distribution to the Fisher information is approximately

$\displaystyle\text{\boldmath$F$}(\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])$ $\displaystyle\approx$
$\displaystyle-\frac{\partial^{2}\log(f_{Y}(y\mid\text{\boldmath$\vartheta$}))}{\partial\text{\boldmath$\vartheta$}\partial\text{\boldmath$\vartheta$}^{\top}}$
$\displaystyle=$
$\displaystyle-\left(\text{\boldmath$a$}(y)\text{\boldmath$a$}(y)^{\top}\left\\{\frac{f_{Z}^{\prime\prime}(\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$})}{f_{Z}(\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$})}-\left[\frac{f_{Z}^{\prime}(\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$})}{f_{Z}(\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$})}\right]^{2}\right\\}\right.$
$\displaystyle\left.\quad-\frac{\text{\boldmath$a$}^{\prime}(y){\text{\boldmath$a$}^{\prime}(y)}^{\top}}{{(\text{\boldmath$a$}^{\prime}(y)}^{\top}\text{\boldmath$\vartheta$})^{2}}\right)$

(NB: the weight to $\text{\boldmath$a$}(y)\text{\boldmath$a$}(y)^{\top}$ is
constant one for $F_{Z}=\Phi$). For a censored or discrete observation, we
have the following contribution to the Fisher information

$\displaystyle\text{\boldmath$F$}(\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])$ $\displaystyle=$
$\displaystyle-\frac{\partial^{2}\log(\mathcal{L}(\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}]))}{\partial\text{\boldmath$\vartheta$}\partial\text{\boldmath$\vartheta$}^{\top}}$
$\displaystyle=$
$\displaystyle-\left\\{\frac{f_{Z}^{\prime}(\text{\boldmath$a$}(\bar{y})^{\top}\text{\boldmath$\vartheta$})\text{\boldmath$a$}(\bar{y})\text{\boldmath$a$}(\bar{y})^{\top}-f_{Z}^{\prime}(\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}\text{\boldmath$\vartheta$})\text{\boldmath$a$}(\underaccent{\bar}{\ry})\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}}{F_{Z}(\text{\boldmath$a$}(\bar{y})^{\top}\text{\boldmath$\vartheta$})-F_{Z}(\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}\text{\boldmath$\vartheta$})}\right.$
$\displaystyle\quad-\frac{[f_{Z}(\text{\boldmath$a$}(\bar{y})^{\top}\text{\boldmath$\vartheta$})\text{\boldmath$a$}(\bar{y})-f_{Z}(\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}\text{\boldmath$\vartheta$})\text{\boldmath$a$}(\underaccent{\bar}{\ry})]}{[F_{Z}(\text{\boldmath$a$}(\bar{y})^{\top}\text{\boldmath$\vartheta$})-F_{Z}(\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}\text{\boldmath$\vartheta$}]^{2}}\times$
$\displaystyle\left.\qquad[f_{Z}(\text{\boldmath$a$}(\bar{y})^{\top}\text{\boldmath$\vartheta$})\text{\boldmath$a$}(\bar{y})^{\top}-f_{Z}(\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}\text{\boldmath$\vartheta$})\text{\boldmath$a$}(\underaccent{\bar}{\ry})^{\top}]\right\\}.$

For a truncated observation, the Fisher information is given by
$\text{\boldmath$F$}(\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])-\text{\boldmath$F$}(\text{\boldmath$\vartheta$}\mid
Y\in(y_{l},y_{r}])$.

We will first discuss the asymptotic properties of the maximum likelihood
estimator $\hat{\text{\boldmath$\vartheta$}}_{N}$ in the parametric setting
with fixed parameters $\vartheta$ in both the discrete and continuous case.
For continuous variables $Y$ and a transformation function parameterised using
a Bernstein polynomial, results for sieve maximum likelihood estimation, where
the number of parameters increases with $N$, are then discussed in Subsection
3.2.

###  3.1 Parametric Inference

Conditions on the densities of the error distribution $f_{Z}$ and the basis
functions $a$ ensuring consistency and asymptotic normality of the sequence of
maximum likelihood estimators $\hat{\text{\boldmath$\vartheta$}}_{N}$ and an
estimator of their asymptotic covariance matrix are given in the following
three theorems. Due to the full parameterisation of the model, the proofs are
simple standard results for likelihood asymptotics, and a more complex
analysis (as required for estimation equations in the presence of a nuisance
parameter $h_{Y}$, for example in Cheng _et al._ , 1995; Chen _et al._ , 2002)
is not necessary. We will restrict ourselves to absolutely continuous or
discrete random variables $Y$, where the likelihood is given in terms of the
density $f_{Y}(y\mid\text{\boldmath$\vartheta$})$. Furthermore, we will only
study the case of a correctly specified transformation
$h=\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$ and refer the reader
to Hothorn _et al._ (2014), where consistency results for arbitrary $h$ are
given.

######  Theorem 1.

For
$Y_{1},\dots,Y_{N}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\mathbb{P}_{Y,\text{\boldmath$\vartheta$}_{0}}$
and under the assumptions (A1) the parameter space $\Theta$ is compact and
(A2)
$\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}[\sup_{\text{\boldmath$\vartheta$}\in\Theta}|\log(f_{Y}(Y\mid\text{\boldmath$\vartheta$}))|]<\infty$
where $\text{\boldmath$\vartheta$}_{0}$ is well-separated:

$\displaystyle\sup_{\text{\boldmath$\vartheta$};|\text{\boldmath$\vartheta$}-\text{\boldmath$\vartheta$}_{0}|\geq\epsilon}\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}[\log(f_{Y}(Y\mid\text{\boldmath$\vartheta$}))]<\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}[\log(f_{Y}(Y\mid\text{\boldmath$\vartheta$}_{0}))],$

the sequence of estimators $\hat{\text{\boldmath$\vartheta$}}_{N}$ converges
to $\text{\boldmath$\vartheta$}_{0}$ in probability,
$\hat{\text{\boldmath$\vartheta$}}_{N}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\text{\boldmath$\vartheta$}_{0}$,
as $N\to\infty$.

###### Proof.

The log-likelihood is continuous in $\vartheta$, and due to (A2), each log-
likelihood contribution is dominated by an integrable function. Thus, the
result follows from van der Vaart (1998) (Theorem 5.8 with Example 19.7; see
note at bottom of page 46). ∎

######  Remark 1.

Assumption (A1) is made for convenience, and relaxations of such a condition
are given in van de Geer (2000) or van der Vaart (1998). The assumptions in
(A2) are rather weak: the first one holds if the functions $a$ are not
arbitrarily ill-posed, and the second one holds if the function
$\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}[\log(f_{Y}(Y\mid\text{\boldmath$\vartheta$}))]$
is strictly convex in $\vartheta$ (if the assumption would not hold, we would
still have convergence to the set
$\mbox{argmax}_{\text{\boldmath$\vartheta$}}\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}[\log(f_{Y}(Y\mid\text{\boldmath$\vartheta$}))]$).

######  Theorem 2.

Under the assumptions of Theorem 1 and in addition (A3)

$\displaystyle\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}\left(\sup_{\text{\boldmath$\vartheta$}}\left|\left|\frac{\partial\log
f_{Y}(Y\mid\text{\boldmath$\vartheta$})}{\partial\text{\boldmath$\vartheta$}}\right|\right|^{2}\right)<\infty,$

(A4)
$\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}(\text{\boldmath$a$}(Y)\text{\boldmath$a$}(Y)^{\top})$
and (for the absolutely continuous case $\mu=\mu_{L}$ only)
$\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}(\text{\boldmath$a$}^{\prime}(Y)\text{\boldmath$a$}^{\prime}(Y)^{\top})$
are nonsingular, and (A5) $0<f_{Z}<\infty$, $\sup|f_{Z}^{\prime}|<\infty$ and
$\sup|f_{Z}^{\prime\prime}|<\infty$, the sequence
$\sqrt{N}(\hat{\text{\boldmath$\vartheta$}}_{N}-\text{\boldmath$\vartheta$}_{0})$
is asymptotically normal with mean zero and covariance matrix

$\displaystyle\Sigma_{\text{\boldmath$\vartheta$}_{0}}=\left(\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}\left(-\frac{\partial^{2}\log
f_{Y}(Y\mid\text{\boldmath$\vartheta$})}{\partial\text{\boldmath$\vartheta$}\partial\text{\boldmath$\vartheta$}^{\top}}\right)\right)^{-1},$

as $N\to\infty$.

###### Proof.

Because the map
$\text{\boldmath$\vartheta$}\mapsto\sqrt{f_{Y}(y\mid\text{\boldmath$\vartheta$})}$
is continuously differentiable in $\vartheta$ for all $y$ in both the discrete
and absolutely continuous case and the matrix

$\displaystyle\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}\left(\left[\frac{\partial\log
f_{Y}(Y\mid\text{\boldmath$\vartheta$})}{\partial\text{\boldmath$\vartheta$}}\right]\left[\frac{\partial\log
f_{Y}(Y\mid\text{\boldmath$\vartheta$})}{\partial\text{\boldmath$\vartheta$}}\right]^{\top}\right)$

is continuous in $\vartheta$ as given in (4) and (5), the transformation
family $\mathbb{P}_{Y,\Theta}$ is differentiable in quadratic mean with Lemma
7.6 in van der Vaart (1998). Furthermore, assumptions (A4-5) ensure that the
expected Fisher information matrix is nonsingular at
$\text{\boldmath$\vartheta$}_{0}$. With the consistency and (A3), the result
follows from Theorem 5.39 in van der Vaart (1998). ∎

######  Remark 2.

Assumption (A4) is valid for the densities $f_{Z}$ of the normal, logistic and
minimum extreme value distribution. The Fisher information (3) and (3)
evaluated at the maximum likelihood estimator
$\hat{\text{\boldmath$\vartheta$}}_{N}$ can be used to estimate the covariance
matrix $\Sigma_{\text{\boldmath$\vartheta$}_{0}}$.

######  Theorem 3.

Under the assumptions of Theorem 2 and assuming
$\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}|\text{\boldmath$F$}(\text{\boldmath$\vartheta$}_{0}\mid
Y)|<\infty$, a consistent estimator for
$\Sigma_{\text{\boldmath$\vartheta$}_{0}}$ is given by

$\displaystyle\hat{\Sigma}_{\text{\boldmath$\vartheta$}_{0},N}=\left(N^{-1}\sum_{i=1}^{N}\text{\boldmath$F$}(\hat{\text{\boldmath$\vartheta$}}_{N}\mid
Y_{i})\right)^{-1}.$

###### Proof.

With the law of large numbers we have

$\displaystyle
N^{-1}\sum_{i=1}^{N}\text{\boldmath$F$}(\text{\boldmath$\vartheta$}_{0}\mid
Y_{i})=N^{-1}\sum_{i=1}^{N}-\frac{\partial^{2}\log
f_{Y}(Y_{i}\mid\text{\boldmath$\vartheta$})}{\partial\text{\boldmath$\vartheta$}\partial\text{\boldmath$\vartheta$}^{\top}}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}\left(-\frac{\partial^{2}\log
f_{Y}(Y\mid\text{\boldmath$\vartheta$})}{\partial\text{\boldmath$\vartheta$}\partial\text{\boldmath$\vartheta$}^{\top}}\right)=\Sigma_{\text{\boldmath$\vartheta$}_{0}}^{-1}.$

Because the map
$\text{\boldmath$\vartheta$}\mapsto\text{\boldmath$F$}(\text{\boldmath$\vartheta$}\mid
y)$ is continuous for all $y$ (as can be seen from (3) and (3)), the result
follows with Theorem 1. ∎

Based on Theorems 1-3, we can perform standard likelihood inference on the
model parameters $\vartheta$. In particular, we can construct confidence
intervals and confidence bands for the conditional distribution function from
confidence intervals and bands for the linear functions
$\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$. We complete this part
by formally defining the class of transformation models.

######  Definition 4 (Transformation model).

The triple $(F_{Z},\text{\boldmath$a$},\text{\boldmath$\vartheta$})$ is called
transformation model.

The transformation model
$(F_{Z},\text{\boldmath$a$},\text{\boldmath$\vartheta$})$ fully defines the
distribution of $Y$ via
$F_{Y}=F_{Z}\circ\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$ and
thus the corresponding likelihood
$\mathcal{L}(\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])$. Our definition of transformation
models as $(F_{Z},\text{\boldmath$a$},\text{\boldmath$\vartheta$})$ is
strongly tied to the idea of structural inference (Fraser, 1968) and group
families (Lehmann, 1983) or group models (Bickel _et al._ , 1993). Fraser
(1968) described a measurement model $\mathbb{P}_{Y}$ for $Y$ by an error
distribution $\mathbb{P}_{Z}$ and a structural equation $Y=g\circ Z$, where
$g$ is a linear function, thereby extending the location-scale family
$Y=\alpha+\sigma Z$ introduced by Fisher (1934) and refined by Pitman (1939).
Group models consist of distributions generated by possibly non-linear $g$.
The main difference to these classical approaches is that we parameterise $h$
instead of $g=h^{-1}$. By extending the linear transformation functions $g$
dealt with by Fraser (1968) to non-linear transformations, we approximate the
potentially non-linear transformation functions $h=g^{-1}=F_{Z}^{-1}\circ
F_{Y}$ by $\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$, with
subsequent estimation of the parameters $\vartheta$. For given parameters
$\vartheta$, a sample from $\mathbb{P}_{Y}$ can be drawn by the probability
integral transform, i.e.
$Z_{1},\dots,Z_{N}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\mathbb{P}_{Z}$
is drawn and then
$Y_{i}=\inf\\{y\in\Xi\mid\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}\geq
Z_{i}\\}$. This generalises the method published by Bender _et al._ (2005)
from the Cox model to all conditional transformation models.

###  3.2 Non-parametric Inference

For continuous responses $Y$, any unknown transformation $h$ can be
approximated by Bernstein polynomials of increasing order (Farouki, 2012). For
uncensored and right-censored responses and under the same conditions for
$F_{Z}$ as stated in Subsection 3.1, McLain and Ghosh (2013) showed that the
non-parametric sieve maximum likelihood estimator is consistent with rate of
convergence $N^{2/5}$ for $h$ with continuous bounded second derivatives in
unconditional and linear transformation models (see Subsection 4.3). In the
latter class, the linear shift parameters $\beta$ are asympotically normal and
semi-parametrically efficient. Numerical approximations to the observed Fisher
information $\text{\boldmath$F$}(\hat{\text{\boldmath$\vartheta$}}_{N}\mid
Y\in(\underaccent{\bar}{\ry},\bar{y}])$ were shown to lead to appropriate
standard errors of $\hat{\text{\boldmath$\beta$}}_{N}$ by McLain and Ghosh
(2013). Hothorn _et al._ (2014) established the consistency of boosted non-
parametric conditional transformation models (see Subsection 4.2). For sieve
maximum likelihood estimation in the class of conditional transformation
models, the techniques employed by McLain and Ghosh (2013) require minor
technical extensions, which are omitted here.

In summary, the same limiting distribution arises under both the parametric
and the non-parametric paradigm for transformation functions parameterised or
approximated using Bernstein polynomials, respectively. In the latter case,
the target is then the best approximated transformation function with
Bernstein polynomials, say $h_{N}^{\star}$ (where the index $N$ indicates that
we use a more complex approximation when $N$ increases). If the approximation
error $h_{N}^{\star}-h$ is of smaller order than the convergence rate of the
estimator, the estimator’s target becomes the true underlying transformation
function $h$, and otherwise a bias for estimating $h$ remains.

##  4 Applications

The definition of transformation models tailored for specific situations
“only” requires the definition of a suitable basis function $a$ and a choice
of $F_{Z}$. In this section, we will discuss specific transformation models
for unconditional and conditional distributions of ordered and unordered
categorical, discrete and continuous responses $Y$. Note that the likelihood
function $\mathcal{L}$ allows all these models to be fitted to arbitrarily
censored or truncated responses; for brevity, we will not elaborate on the
details.

###  4.1 Unconditional Transformation Models

##### Finite Sample Space

For ordered categorical responses $Y$ from a finite sample space
$\Xi=\\{y_{1},\dots,y_{K}\\}$, we assign one parameter to each element of the
sample space except $y_{K}$. This corresponds to the basis function
$\text{\boldmath$a$}(y_{k})=\text{\boldmath$e$}_{K-1}(k)$, where
$\text{\boldmath$e$}_{K-1}(k)$ is the unit vector of length $K-1$, with its
$k$th element being one. The transformation function $h$ is

$\displaystyle h(y_{k})$ $\displaystyle=$
$\displaystyle\text{\boldmath$e$}_{K-1}(k)^{\top}\text{\boldmath$\vartheta$}=\vartheta_{k}\in\mathbb{R},\quad
1\leq k<K,\quad\text{st}\quad\vartheta_{1}<\dots<\vartheta_{K\textendash 1},$

with $h(y_{K})=\infty$, and the unconditional distribution function of $F_{Y}$
is $F_{Y}(y_{k})=F_{Z}(\vartheta_{k})$. This parameterisation underlies the
common proportional odds and proportional hazards model for ordered
categorical data (Tutz, 2012). Note that monotonicity of $h$ is guaranteed by
the $K-2$ linear constraints
$\vartheta_{2}-\vartheta_{1}>0,\dots,\vartheta_{K-1}-\vartheta_{K-2}>0$ when
constrained optimisation is performed. In the absence of censoring or
truncation and with $\vartheta_{0}=-\infty,\vartheta_{K}=\infty$, we obtain
the maximum likelihood estimator for $\vartheta$ as

$\displaystyle\hat{\text{\boldmath$\vartheta$}}_{N}$ $\displaystyle=$
$\displaystyle\operatorname*{{arg\,max}}_{\vartheta_{1}<\dots<\vartheta_{K-1}}\sum_{i=1}^{N}\log(F_{Z}(\vartheta_{k(i)})-F_{Z}(\vartheta_{k(i)-1}))$
$\displaystyle=$
$\displaystyle\left(F_{Z}^{-1}\left(N^{-1}\sum_{i=1}^{N}\mathds{1}(Y_{i}\leq
y_{1})\right),\dots,F_{Z}^{-1}\left(N^{-1}\sum_{i=1}^{N}\mathds{1}(Y_{i}\leq
y_{K-1})\right)\right)^{\top}$

because $\hat{\pi}_{k}=N^{-1}\sum_{i=1}^{N}\mathds{1}(Y_{i}=y_{k}),1\leq k<K$
maximises the equivalent multinomial (or empirical) log-likelihood
$\sum_{i=1}^{N}\log(\pi_{k(i)})$, and we can rewrite this estimator as

$\displaystyle\hat{\pi}_{k}=N^{-1}\left(\sum_{i=1}^{N}\mathds{1}(Y_{i}\leq
y_{k})-\mathds{1}(k>1)\sum_{i=1}^{N}\mathds{1}(Y_{i}\leq y_{k-1})\right),1\leq
k<K.$

The estimated distribution function $\hat{F}_{Y,N}=F_{Z}\circ\hat{h}_{N}$ is
invariant with respect to $F_{Z}$. Assumption (A4) is valid for these basis
functions because we have
$\mathbb{E}_{\text{\boldmath$\vartheta$}_{0}}(\text{\boldmath$e$}_{K-1}(Y)\text{\boldmath$e$}_{K-1}(Y)^{\top})=\operatorname{diag}(\mathbb{P}(Y=y_{k})),1\leq
k<K$ for $Y\sim\mathbb{P}_{Y,\text{\boldmath$\vartheta$}_{0}}$.

If we define the sample space $\Xi$ as the set of unique observed values and
the probability measure as the empirical cumulative distribution function
(ECDF), putting mass $N^{-1}$ on each observation, we see that this particular
parameterisation is equivalent to an empirical likelihood approach and we get
$\hat{h}_{N}=F_{Z}^{-1}\circ\text{ECDF}$. Note that although the
transformation function depends on the choice of $F_{Z}$, the estimated
distribution function $\hat{F}_{Y,N}=F_{Z}\circ\hat{h}_{N}=\text{ECDF}$ does
not and is simply the non-parametric empirical maximum likelihood estimator. A
smoothed version of this estimator for continuous responses is discussed in
the next paragraph.

##### Infinite Sample Space

For continuous responses $Y$, the parameterisation
$h(y)=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}$, and thus also
$\hat{F}_{Y,N}$, should be smooth in $y$; therefore, any polynomial or spline
basis is a suitable choice for $a$. For the empirical experiments in Section
5, we applied Bernstein polynomials (for an Overview, see Farouki, 2012) of
order $M$ ($P=M+1$) defined on the interval
$[\underaccent{\bar}{\imath},\bar{\imath}]$ with

$\displaystyle\text{\boldmath$a$}_{\text{Bs},M}(y)$ $\displaystyle=$
$\displaystyle(M+1)^{-1}(f_{\text{Be}(1,M+1)}(\tilde{y}),\dots,f_{\text{Be}(m,M-m+1)}(\tilde{y}),\dots,f_{\text{Be}(M+1,1)}(\tilde{y}))^{\top}\in\mathbb{R}^{M+1}$
$\displaystyle h(y)$ $\displaystyle=$
$\displaystyle\text{\boldmath$a$}_{\text{Bs},M}(y)^{\top}\text{\boldmath$\vartheta$}=\sum_{m=0}^{M}\vartheta_{m}f_{\text{Be}(m+1,M-m+1)}(\tilde{y})/(M+1)$
$\displaystyle h^{\prime}(y)$ $\displaystyle=$
$\displaystyle\text{\boldmath$a$}_{\text{Bs},M}^{\prime}(y)^{\top}\text{\boldmath$\vartheta$}=\sum_{m=0}^{M-1}(\vartheta_{m+1}-\vartheta_{m})f_{\text{Be}(m+1,M-m)}(\tilde{y})M/((M+1)(\bar{\imath}-\underaccent{\bar}{\imath})),$

where
$\tilde{y}=(y-\underaccent{\bar}{\imath})/(\bar{\imath}-\underaccent{\bar}{\imath})\in[0,1]$
and $f_{\text{Be}(m,M)}$ is the density of the Beta distribution with
parameters $m$ and $M$. This choice is computationally attractive because
strict monotonicity can be formulated as a set of $M$ linear constraints on
the parameters $\vartheta_{m}<\vartheta_{m+1}$ for all $m=0,\dots,M$ (Curtis
and Ghosh, 2011). Therefore, application of constrained optimisation
guarantees monotonic estimates $\hat{h}_{N}$. The basis contains an intercept.
We obtain smooth plug-in estimators for the distribution, density, hazard and
cumulative hazard functions as
$\hat{F}_{Y,N}=F_{Z}\circ\text{\boldmath$a$}_{\text{Bs},M}^{\top}\hat{\text{\boldmath$\vartheta$}}_{N}$,
$\hat{f}_{Y,N}=f_{Z}\circ\text{\boldmath$a$}_{\text{Bs},M}^{\top}\hat{\text{\boldmath$\vartheta$}}_{N}\times{\text{\boldmath$a$}_{\text{Bs},M}^{\prime}}^{\top}\hat{\text{\boldmath$\vartheta$}}_{N}$,
$\hat{\lambda}_{Y,N}=\hat{f}_{Y,N}/(1-\hat{F}_{Y,N})$ and
$\hat{\Lambda}_{Y,N}=-\log(1-\hat{F}_{Y,N})$. The estimator
$\hat{F}_{Y,N}=F_{Z}\circ\text{\boldmath$a$}_{\text{Bs},M}^{\top}\hat{\text{\boldmath$\vartheta$}}_{N}$
must not be confused with the estimator
$\hat{F}_{Y,N}=\text{\boldmath$a$}_{\text{Bs},M}^{\top}\hat{\text{\boldmath$p$}}$
for $Y\in[0,1]$ obtained from the smoothed empirical distribution function
with coefficients
$\hat{\text{\boldmath$p$}}_{m+1}=\sum_{i=1}^{N}\mathds{1}(Y_{i}\leq m/M)/N$
corresponding to probabilities evaluated at the quantiles $m/M$ for
$m=0,\dots,M$ (Babu _et al._ , 2002).

The question arises how the degree of the polynomial affects the estimated
distribution function. On the one hand, the model
$(\Phi,\text{\boldmath$a$}_{\text{Bs},1},\text{\boldmath$\vartheta$})$ only
allows linear transformation functions of a standard normal and $F_{Y}$ is
restricted to the normal family. On the other hand,
$(\Phi,\text{\boldmath$a$}_{\text{Bs},N-1},\text{\boldmath$\vartheta$})$ has
one parameter for each observation and $\hat{F}_{Y,N}$ is the non-parametric
maximum likelihood estimator ECDF, which, by the Glivenko-Cantelli lemma,
converges to $F_{Y}$. In this sense, we cannot choose a “too large” value for
$M$. This is a consequence of the monotonicity constraint on the estimator
$\text{\boldmath$a$}^{\top}\hat{\text{\boldmath$\vartheta$}}_{N}$, which, in
this extreme case, just interpolates the step function
$F_{Z}^{-1}\circ\text{ECDF}$. Empirical evidence for the insensitivity of
results when $M$ is large can be found in Hothorn (2017b) and in the
discussion.

###  4.2 Conditional Transformation Models

In the following, we will discuss a cascade of increasingly complex
transformation models where the transformation function $h$ may depend on
explanatory variables $\text{\boldmath$X$}\in\chi$. We are interested in
estimating the conditional distribution of $Y$ given
$\text{\boldmath$X$}=\text{\boldmath$x$}$. The corresponding distribution
function $F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}$ can be written as
$F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=F_{Z}(h(y\mid\text{\boldmath$x$}))$.
The transformation function
$h(\cdot\mid\text{\boldmath$x$}):\Xi\rightarrow\mathbb{R}$ is said to be
conditional on $x$. Following the arguments presented in the proof of
Corollary 1, it is easy to see that for each $x$, there exists a strictly
monotonic transformation function
$h(\cdot\mid\text{\boldmath$x$})=F_{Z}^{-1}\circ
F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}$ such that
$F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=F_{Z}(h(y\mid\text{\boldmath$x$}))$.
Because this class of conditional transformation models and suitable
parameterisations were introduced by Hothorn _et al._ (2014), we will only
sketch the most important aspects here.

Let $\text{\boldmath$b$}:\chi\rightarrow\mathbb{R}^{Q}$ denote a basis
transformation of the explanatory variables. The joint basis for both $y$ and
$x$ is called
$\text{\boldmath$c$}:\Xi\times\chi\rightarrow\mathbb{R}^{d(P,Q)}$; its
dimension $d(P,Q)$ depends on the way the two basis functions $a$ and $b$ are
combined (e.g.
$\text{\boldmath$c$}=(\text{\boldmath$a$}^{\top},\text{\boldmath$b$}^{\top})^{\top}\in\mathbb{R}^{P+Q}$
or
$\text{\boldmath$c$}=(\text{\boldmath$a$}^{\top}\otimes\text{\boldmath$b$}^{\top})^{\top}\in\mathbb{R}^{PQ}$).
The conditional transformation function is now parameterised as
$h(y\mid\text{\boldmath$x$})=\text{\boldmath$c$}(y,\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$}$.
One important special case is the simple transformation function
$h(y\mid\text{\boldmath$x$})=h_{Y}(y)+h_{\text{\boldmath$x$}}(\text{\boldmath$x$})$,
where the explanatory variables only contribute a shift
$h_{\text{\boldmath$x$}}(\text{\boldmath$x$})$ to the conditional
transformation function. Often this shift is assumed to be linear in $x$;
therefore, we use the function
$m(\text{\boldmath$x$})=\text{\boldmath$b$}(\text{\boldmath$x$})^{\top}\text{\boldmath$\beta$}=\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$
to denote linear shifts. Here
$\text{\boldmath$b$}(\text{\boldmath$x$})=\tilde{\text{\boldmath$x$}}$ is one
row of the design matrix without intercept. These simple models correspond to
the joint basis
$\text{\boldmath$c$}(y,\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$}=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{1}+\text{\boldmath$b$}(\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$}_{2}$,
with $h_{Y}(y)=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{1}$
and
$h_{\text{\boldmath$x$}}(\text{\boldmath$x$})=\text{\boldmath$b$}(\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$}_{2}=m(\text{\boldmath$x$})=\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$.
The results presented in Section 3, including Theorems 1, 2 and 3, carry over
in the fixed design case when $a$ is replaced by $c$.

In the rest of this section, we will present classical models that can be
embedded in the larger class of conditional transformation models and some
novel models that can be implemented in this general framework.

###  4.3 Classical Transformation Models

##### Linear Model

The normal linear regression model
$Y\sim\operatorname{N}(\alpha+m(\text{\boldmath$x$}),\sigma^{2})$ with
conditional distribution function

$\displaystyle
F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=\Phi\left(\frac{y-(\alpha+m(\text{\boldmath$x$}))}{\sigma}\right)$

can be understood as a transformation model with transformation function
$h(y\mid\text{\boldmath$x$})=y/\sigma-\alpha/\sigma-m(\text{\boldmath$x$})/\sigma$
parameterised via basis functions
$\text{\boldmath$a$}(y)=(y,1)^{\top},\text{\boldmath$b$}(\text{\boldmath$x$})=\tilde{\text{\boldmath$x$}}$
and
$\text{\boldmath$c$}=(\text{\boldmath$a$}^{\top},\text{\boldmath$b$}^{\top})^{\top}$
with parameters
$\text{\boldmath$\vartheta$}=(\sigma^{-1},-\sigma^{-1}\alpha,-\sigma^{-1}\text{\boldmath$\beta$}^{\top})^{\top}$
under the constraint $\sigma>0$, or in more compact notation
$(\Phi,(y,1,\tilde{\text{\boldmath$x$}}^{\top})^{\top},\text{\boldmath$\vartheta$})$.
The parameters of the model are the inverse standard deviation and the inverse
negative coefficient of variation instead of the mean and variance of the
original normal distribution. For “exact continuous” observations, the
likelihood $\mathcal{L}$ is equivalent to least-squares, which can be
maximised with respect to $\alpha$ and $\beta$ without taking $\sigma$ into
account. This is not possible for censored or truncated observations, where we
need to evaluate the conditional distribution function that depends on all
parameters; this model is called Type I Tobit model (Tobin, 1958) (although
only the likelihood changes under censoring and truncation, but the model does
not). Using an alternative basis function $c$ would allow arbitrary non-normal
conditional distributions of $Y$ and the simple shift model
$\text{\boldmath$c$}(y,\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$}=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{1}+\text{\boldmath$b$}(\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$}_{2}$
is then a generalisation of additive models and leads to the interpretation

$\displaystyle\mathbb{E}_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(\text{\boldmath$a$}(Y)^{\top}\text{\boldmath$\vartheta$}_{1})=-\text{\boldmath$b$}(\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$}_{2}.$

The choice $\text{\boldmath$a$}=(1,\log)^{\top}$ implements the log-normal
model for $Y>0$. Implementation of a Bernstein basis
$\text{\boldmath$a$}=\text{\boldmath$a$}_{\text{Bs},M}$ allows arbitrarily
shaped distributions, i.e. a transition from the normal family to the
transformation family, and thus likelihood inference on
$\text{\boldmath$\vartheta$}_{2}$ without strict assumptions on the
distribution of $Y$. The transformation
$\text{\boldmath$a$}_{\text{Bs},M}(y)^{\top}\text{\boldmath$\vartheta$}_{1}$
must increase monotonically in $y$. Maximisation of the log-likelihood under
the linear inequality constraint
$\text{\boldmath$D$}_{M+1}\text{\boldmath$\vartheta$}_{1}>0$, with
$\text{\boldmath$D$}_{M+1}$ representing first-order differences, implements
this requirement.

##### Continuous “Survival Time” Models

For a continuous response $Y>0$, the model

$\displaystyle
F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=F_{Z}\left(\frac{\log(y)-(\alpha+m(\text{\boldmath$x$}))}{\sigma}\right)$

with basis functions $\text{\boldmath$a$}(y)=(1,\log(y))$ and
$\text{\boldmath$b$}(\text{\boldmath$x$})=\tilde{\text{\boldmath$x$}}$ and
parameters
$\text{\boldmath$\vartheta$}=(-\alpha,\sigma^{-1},-\text{\boldmath$\beta$}^{\top})^{\top}$
under the constraint $\sigma>0$ is called the accelerated failure time (AFT)
model. The model
$(F_{\operatorname{MEV}},(1,\log,\tilde{\text{\boldmath$x$}}^{\top})^{\top},(-\vartheta_{1},1,-\text{\boldmath$\beta$}^{\top})^{\top})$
with $\sigma\equiv 1$ (and thus fixed transformation function $\log$) is the
exponential AFT model because it implies an exponential distribution of $Y$.
When the parameter $\sigma>0$ is estimated from the data, the model
$(F_{\operatorname{MEV}},(1,\log,\tilde{\text{\boldmath$x$}}^{\top})^{\top},\text{\boldmath$\vartheta$})$
is called the Weibull model,
$(F_{\operatorname{SL}},(1,\log,\tilde{\text{\boldmath$x$}}^{\top})^{\top},\text{\boldmath$\vartheta$})$
is the log-logistic AFT model and
$(\Phi,(1,\log,\tilde{\text{\boldmath$x$}}^{\top})^{\top},\text{\boldmath$\vartheta$})$
is the log-normal AFT model. For a continuous (not necessarily positive)
response $Y$, the model
$F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=F_{\operatorname{MEV}}(h_{Y}(y)-m(\text{\boldmath$x$}))$
is called the proportional hazards, relative risk or Cox model. The
transformation function $h_{Y}$ equals the log-cumulative baseline hazard and
is treated as a nuisance parameter in the partial likelihood framework, where
only the regression coefficients $\beta$ are estimated. Given
$\hat{\text{\boldmath$\beta$}}$, non-parametric maximum likelihood estimators
are typically applied to obtain $\hat{h}_{Y}$. Here, we parameterise this
function as
$h_{Y}(y)=\log(\Lambda_{Y}(y))=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{1}$
(for example, using $\text{\boldmath$a$}=\text{\boldmath$a$}_{\text{Bs},M}$)
and fit all parameters in the model
$(F_{\operatorname{MEV}},(\text{\boldmath$a$}^{\top},\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\text{\boldmath$\vartheta$}_{1}^{\top},-\text{\boldmath$\beta$}^{\top})^{\top})$
simultaneously. The model is highly popular because $m(\text{\boldmath$x$})$
is the log-hazard ratio to $m(\mathbf{0})$. For the special case of right-
censored survival times, this parameterisation of the Cox model was studied
theoretically and empirically by McLain and Ghosh (2013). Changing the
distribution function in the Cox model from $F_{\operatorname{MEV}}$ to
$F_{\operatorname{SL}}$ results in the proportional odds model
$(F_{\operatorname{SL}},(\text{\boldmath$a$}^{\top},\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\text{\boldmath$\vartheta$}_{1}^{\top},-\text{\boldmath$\beta$}^{\top})^{\top})$;
its name comes from the interpretation of $m(\text{\boldmath$x$})$ as the
constant log-odds ratio of the odds
$O_{Y}(y\mid\text{\boldmath$X$}=\text{\boldmath$x$})$ and
$O_{Y}(y\mid\text{\boldmath$x$}=\mathbf{0})$. An additive hazards model with
the conditional hazard function
$\lambda_{Y}(y\mid\text{\boldmath$X$}=\text{\boldmath$x$})=\lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$
results from the choice $F_{Z}(z)=F_{\operatorname{Exp}}(z)=1-\exp(-z)$
(Aranda-Ordaz, 1983) under the additional constraint
$\lambda_{Y}(y\mid\text{\boldmath$X$}=\text{\boldmath$x$})>0$. In this case,
the function $\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{1}>0$
is the positive baseline cumulative hazard function
$\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$.

##### Discrete Models

For ordered categorical responses $y_{1}<\dots<y_{K}$, the conditional
distribution
$F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y_{k})=F_{Z}(h_{Y}(y_{k})-m(\text{\boldmath$x$}))$
is a transformation model with
$\text{\boldmath$a$}(y_{k})=\text{\boldmath$e$}_{K-1}(k)$. The model
$(F_{\operatorname{SL}},(\text{\boldmath$a$}^{\top},\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\text{\boldmath$\vartheta$}_{1}^{\top},-\text{\boldmath$\beta$}^{\top})^{\top})$
is called the discrete proportional odds model and
$(F_{\operatorname{MEV}},(\text{\boldmath$a$}^{\top},\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\text{\boldmath$\vartheta$}_{1}^{\top},-\text{\boldmath$\beta$}^{\top})^{\top})$
is the discrete proportional hazards model. Here, $m(\text{\boldmath$x$})$ is
the log-odds ratio or log-hazard ratio to $m(\mathbf{0})$ independent of $k$;
details are given in Tutz (2012). For the special case of a binary response
($K=2$), the transformation model
$(F_{\operatorname{SL}},(\mathds{1}(k=1),\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\vartheta_{1},-\text{\boldmath$\beta$}^{\top})^{\top})$
is the logistic regression model,
$(\Phi,(\mathds{1}(k=1),\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\vartheta_{1},-\text{\boldmath$\beta$}^{\top})^{\top})$
is the probit model and
$(F_{\operatorname{MEV}},(\mathds{1}(k=1),\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\vartheta_{1},-\text{\boldmath$\beta$}^{\top})^{\top})$
is called the complementary log-log model. Note that the transformation
function $h_{Y}$ is given by the basis function
$\text{\boldmath$a$}=\mathds{1}(k=1)$, i.e. $\vartheta_{1}$ is just the
intercept. The connection between standard binary regression models and
transformation models is explained in more detail by Doksum and Gasko (1990).

##### Linear Transformation Model

The transformation model
$(F_{Z},(\text{\boldmath$a$}^{\top},\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\text{\boldmath$\vartheta$}_{1}^{\top},-\text{\boldmath$\beta$}^{\top})^{\top})$
for any $a$ and $F_{Z}$ is called the linear transformation model and contains
all models discussed in this subsection. Note that the transformation of the
response
$h_{Y}(y)=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{1}$ is
non-linear in all models of interest (AFT, Cox, etc.), and the term “linear”
only refers to a linear shift $m(\text{\boldmath$x$})$ of the explanatory
variables. Partially linear or additive transformation models allow non-linear
shifts as part of a partially smooth basis $b$, i.e. in the form of an
additive model. The number of constraints only depends on the basis $a$ but
not on the explanatory variables.

###  4.4 Extension of Classical Transformation Models

A common property of all classical transformation models is the additivity of
the response transformation and the shift, i.e. the decomposition
$h(y\mid\text{\boldmath$x$})=h_{Y}(y)+h_{\text{\boldmath$x$}}(\text{\boldmath$x$})$
of the conditional transformation function. This assumption is relaxed by the
following extensions of the classical models. Allowing for deviations from
this simple model is also the key aspect for the development of novel
transformation models in the rest of this section.

##### Discrete Non-Proportional Odds and Hazards Models

For ordered categorical responses, the model
$F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y_{k})=F_{Z}(h_{Y}(y_{k})-m_{k}(\text{\boldmath$x$}))$
allows a category-specific shift
$m_{k}(\text{\boldmath$x$})=\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}_{k}$;
with $F_{\operatorname{SL}}$, this cumulative model is called the non-
proportional odds model, and with $F_{\operatorname{MEV}}$, it is the non-
proportional hazards model. Both models can be cast into the transformation
model framework by defining the joint basis
$\text{\boldmath$c$}(y_{k},\text{\boldmath$x$})=(\text{\boldmath$a$}(y_{k})^{\top},\text{\boldmath$a$}(y_{k})^{\top}\otimes\text{\boldmath$b$}(\text{\boldmath$x$})^{\top})^{\top}$
as the Kronecker product of the two simple basis functions
$\text{\boldmath$a$}(y_{k})=\text{\boldmath$e$}_{K-1}(k)$ and
$\text{\boldmath$b$}(\text{\boldmath$x$})=\tilde{\text{\boldmath$x$}}$
(assuming that $b$ does not contain an intercept term). Note that the
conditional transformation function $h(y\mid\text{\boldmath$x$})$ includes an
interaction term between $y$ and $x$. It is also worth noting that for
unordered categorical responses $Y\in\Xi=\\{y_{1},\dots,y_{K}\\}$, the
multinomial model can be estimated by the model
$(F_{Z},\text{\boldmath$c$},\text{\boldmath$\vartheta$})$ under any ordering
of the response categories because the corresponding conditional density

$\displaystyle
f_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y_{k})=\left\\{\begin{array}[]{ll}F_{Z}(\text{\boldmath$c$}(y_{k},\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$})&k=1\\\
F_{Z}(\text{\boldmath$c$}(y_{k},\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$})-F_{Z}(\text{\boldmath$c$}(y_{k-1},\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$})&1<k<K\\\
1-F_{Z}(\text{\boldmath$c$}(y_{k-1},\text{\boldmath$x$})^{\top}\text{\boldmath$\vartheta$})&k=K\end{array}\right.$

is invariant with respect to the ordering applied. For
$F_{\operatorname{SL}}$, this model is called the cumulative logit or partial
proportional odds model (Tutz, 2012). The classical multinomial logit (with
$F_{Z}=F_{\operatorname{SL}}$) or probit (with $F_{Z}=\Phi$) models for an
unordered response $Y\in\Xi=\\{y_{1},\dots,y_{K}\\}$ can be written as models
for the density

$\displaystyle\mathbb{P}(Y=y_{1}\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
$\displaystyle=$ $\displaystyle
F_{Z}(\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}_{1})$
$\displaystyle\mathbb{P}(Y=y_{k}\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
$\displaystyle=$ $\displaystyle
F_{Z}(\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}_{k})-F_{Z}(\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}_{k-1})\quad
1<k<K$
$\displaystyle\mathbb{P}(Y=y_{K}\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
$\displaystyle=$ $\displaystyle
1-F_{Z}(\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}_{K\textendash
1}),$

where the parameters $\text{\boldmath$\beta$}_{k}$ correspond to parameters in
a cumulative model $\mathbb{P}(Y\leq
y_{k})=F_{Z}(\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}_{k})$
for any ordering of the sample space. Of course, the parameter estimates
$\hat{\text{\boldmath$\beta$}}_{k}$ change when the ordering changes and
therefore must not be interpreted directly, but the estimated densities
$\hat{\mathbb{P}}(Y=y_{k}\mid\text{\boldmath$X$}=\text{\boldmath$x$})$ are
invariant with respect to the ordering applied.

##### Time-Varying Effects

One often studied extension of the Cox model is
$F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=F_{Z}(h_{Y}(y)-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}(y))$,
where the regression coefficients $\text{\boldmath$\beta$}(y)$ may change with
time $y$. The Cox model is included with
$\text{\boldmath$\beta$}(y)\equiv\text{\boldmath$\beta$}$, and the model is
often applied to check the proportional hazards assumption. With a smooth
parameterisation of time $y$, for example via
$\text{\boldmath$a$}=\text{\boldmath$a$}_{\text{Bs},M}$, and linear basis
$\text{\boldmath$b$}(\text{\boldmath$x$})=\tilde{\text{\boldmath$x$}}$, the
transformation model
$(F_{\operatorname{MEV}},(\text{\boldmath$a$}^{\top},\text{\boldmath$a$}^{\top}\otimes\text{\boldmath$b$}^{\top})^{\top},\text{\boldmath$\vartheta$})$
implements this Cox model with time-varying (linear) effects. This model (with
arbitrary $F_{Z}$) has also been presented in Foresi and Peracchi (1995) and
is called distribution regression in Chernozhukov _et al._ (2013).

###  4.5 Novel Transformation Models

Due to the broadness of the transformation family, it is straightforward to
set up new models for interesting situations by allowing more complex
transformation functions $h(y\mid\text{\boldmath$x$})$. We will illustrate
this possibility for two simple cases – the independent two-sample situation
and regression models for count data. The generic and most complex
transformation model is called the conditional transformation model and is
explained at the end of this section.

##### Beyond Shift Effects

Assume we observe samples from two groups $A$ and $B$ and want to model the
conditional distribution functions $F_{Y\mid\text{\boldmath$X$}=A}(y)$ and
$F_{Y\mid\text{\boldmath$X$}=B}(y)$ of the response $Y$ in the two groups.
Based on this model, it is often interesting to infer whether the two
distributions are equivalent and, if this is not the case, to characterise how
they differ. Using an appropriate basis function $a$ and the basis
$\text{\boldmath$b$}(\text{\boldmath$x$})=(1,\mathds{1}(B))^{\top}$, the model
$(F_{Z},(\text{\boldmath$a$}^{\top}\otimes\text{\boldmath$b$}^{\top})^{\top},\text{\boldmath$\vartheta$})$
parameterises the conditional transformation function as $h(y\mid
A)=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{1}$ and $h(y\mid
B)=h(y\mid
A)+h_{B-A}(y)=\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{1}+\mathds{1}(B)\text{\boldmath$a$}(y)^{\top}\text{\boldmath$\vartheta$}_{2}$.
Clearly, the second term is constant zero ($h_{B-A}(y)\equiv 0$) iff the two
distributions are equivalent
($F_{Y\mid\text{\boldmath$X$}=A}(y)=F_{Y\mid\text{\boldmath$X$}=B}(y)$ for all
$y$). For the deviation function
$h_{B-A}(y)=\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}_{2}$, we can
apply standard likelihood inference procedures for
$\hat{\text{\boldmath$\vartheta$}}_{2}$ to construct a confidence band or use
a test statistic like
$\max(\hat{\text{\boldmath$\vartheta$}}_{2}/\text{se}(\hat{\text{\boldmath$\vartheta$}}_{2}))$
to assess deviations from zero. If there is evidence for a group effect, we
can use the model to check whether the deviation function is constant, i.e.
$h_{B-A}(y)\equiv c\neq 0$. In this case, the simpler model
$(F_{Z},(\text{\boldmath$a$}^{\top},\mathds{1}(B))^{\top},(\text{\boldmath$\vartheta$}_{1}^{\top},-\beta)^{\top})$
with shift $\beta=-\text{\boldmath$\vartheta$}_{2}$ might be easier to
interpret. This model actually corresponds to a normal ANOVA model with
$F_{Z}=\Phi$ and $\text{\boldmath$a$}(y)^{\top}=(1,y)^{\top}$ or the Cox
proportional hazards model with
$(F_{\operatorname{MEV}},(\text{\boldmath$a$}_{\text{Bs},M}^{\top},\mathds{1}(B))^{\top},(\text{\boldmath$\vartheta$}_{1}^{\top},-\beta)^{\top})$.

##### Count Regression “Without Tears”

Simple models for count data $\Xi=\\{0,1,2,\dots\\}$ almost always suffer from
over-dispersion or excess zeros. The linear transformation model

$\displaystyle
F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=F_{Z}(h_{Y}(y)-m(\text{\boldmath$x$}))$

can be implemented using the basis function
$\text{\boldmath$a$}(y)=\text{\boldmath$a$}_{\text{Bs},M}(\lfloor y\rfloor)$,
and then the parameters of the transformation model
$(F_{Z},(\text{\boldmath$a$}^{\top},\tilde{\text{\boldmath$x$}}^{\top})^{\top},\text{\boldmath$\vartheta$})$
are not affected by over- or under-dispersion because higher moments are
handled by $h_{Y}$ independently of the effects of the explanatory variables
$m(\text{\boldmath$x$})$. If there are excess zeros, we can set up a joint
transformation model

$\displaystyle
F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=F_{Z}(h_{Y}(y)-m(\text{\boldmath$x$})+\mathds{1}(y=0)(\alpha_{0}-m_{0}(\text{\boldmath$x$})))$

such that we have a two-components mixture model consisting of the count
distribution

$\displaystyle
F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)=F_{Z}(h_{Y}(y)-m(\text{\boldmath$x$})),y\in\Xi$

and the probability of an excess zero

$\displaystyle
f_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(0)=F_{Z}(h_{Y}(0)-m(\text{\boldmath$x$})+(\alpha_{0}-m_{0}(\text{\boldmath$x$})))=F_{Z}(h_{Y}(0)+\alpha_{0}-\tilde{\text{\boldmath$x$}}^{\top}(\text{\boldmath$\beta$}+\text{\boldmath$\beta$}_{0}))$

when
$m_{0}(\text{\boldmath$x$})=\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}_{0}$.
Hence, the transformation analogue to a hurdle model with hurdle at zero is
the transformation model
$(F_{Z},(\text{\boldmath$a$}^{\top},\tilde{\text{\boldmath$x$}}^{\top},\mathds{1}(y=0),\mathds{1}(y=0)\tilde{\text{\boldmath$x$}}^{\top})^{\top},(\text{\boldmath$\vartheta$}_{1}^{\top},\text{\boldmath$\beta$}^{\top},\alpha_{0},\text{\boldmath$\beta$}_{0}^{\top})^{\top})$.

##### Conditional Transformation Models

When the conditional transformation function is parameterised by multiple
basis functions
$\text{\boldmath$a$}_{j}(y),\text{\boldmath$b$}_{j}(\text{\boldmath$x$}),j=1,\dots
J$ via the joint basis

$\displaystyle\text{\boldmath$c$}=(\text{\boldmath$a$}_{1}^{\top}\otimes\text{\boldmath$b$}_{1}^{\top},\dots,\text{\boldmath$a$}_{J}^{\top}\otimes\text{\boldmath$b$}_{J}^{\top})^{\top},$

models of the class $(\cdot,\text{\boldmath$c$},\text{\boldmath$\vartheta$})$
are called conditional transformation models (CTMs) with $J$ partial
transformation functions parameterised as
$\text{\boldmath$a$}_{j}^{\top}\otimes\text{\boldmath$b$}_{j}^{\top}$ and
include all special cases discussed in this section. It is convenient to
assume monotonicity for each of the partial transformation functions; thus,
the linear constraints for $\text{\boldmath$a$}_{j}$ are repeated for each
basis function in $\text{\boldmath$b$}_{j}$ (detailed descriptions of linear
constraints for different models in this class are available in Hothorn,
2017b). Hothorn _et al._ (2014) introduced this general model class and
proposed a boosting algorithm for the estimation of transformation functions
$h$ for “exact continuous” responses $Y$. In the likelihood framework
presented here, conditional transformation models can be fitted under
arbitrary schemes of censoring and truncation and classical likelihood
inference for the model parameters $\vartheta$ becomes feasible. Of course,
unlike in the boosting context, the number of model terms $J$ and their
complexity are limited in the likelihood world because the likelihood does not
contain any penalty terms that induce smoothness in the $x$-direction. A more
detailed overview on the class of conditional transformation models can be
found in (Möst _et al._ , 2016).

A systematic overview of linear transformation models with potentially
response-varying effects is given in Table 1. Model nomenclature and
interpretation of the corresponding model parameters is mapped to specific
transformation functions $h$ and distribution functions $F_{Z}$. To the best
of our knowledge, models without names have not yet been discussed in the
literature, and their specific properties await closer investigation.

|  |  | $F_{Z}$
---|---|---|---
|  |  | $\Phi$ | $F_{\operatorname{SL}}$ | $F_{\operatorname{Exp}}$ | $F_{\operatorname{MEV}}$
$\Xi$ | $h$ | Meaning of |  |  |  |
$K=2$ | Binary Regression |
| $\mathds{1}(k=1)\vartheta_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$ |  | probit BGLM | logistic BGLM | clog BGLM | cloglog BGLM
|  | $\vartheta_{1}$ |  | $\log(O_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$ | $\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$ | $\log(\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$
|  | $\beta$ |  | log-OR | AH | log-HR
$K>2$ | Polytomous Regression |
| $\text{\boldmath$e$}_{K-1}(k)^{\top}\text{\boldmath$\vartheta$}_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$ |  |  | discrete PO |  | discrete PH
|  | $\text{\boldmath$e$}_{K-1}(k)^{\top}\text{\boldmath$\vartheta$}_{1}$ |  | $\log(O_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$ | $\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$ | $\log(\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$
|  | $\beta$ |  | log-OR | AH | log-HR
| $\text{\boldmath$e$}_{K-1}(k)^{\top}\text{\boldmath$\vartheta$}_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}(k)$ |  |  | non-PO |  | non-PH
|  | $\text{\boldmath$e$}_{K-1}(k)^{\top}\text{\boldmath$\vartheta$}_{1}$ |  | $\log(O_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$ | $\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$ | $\log(\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$
$\mathbb{N}$ | Count Regression |
| $\text{\boldmath$a$}_{\text{Bs},M}(\lfloor y\rfloor)^{\top}\text{\boldmath$\vartheta$}_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$ |  |  |  |  |
|  | $\text{\boldmath$a$}_{\text{Bs},M}(\lfloor y\rfloor)^{\top}\text{\boldmath$\vartheta$}_{1}$ |  | $\log(O_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$ | $\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$ | $\log(\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$
|  | $\beta$ |  | log-OR | AH | log-HR
| $\text{\boldmath$a$}_{\text{Bs},M}(\lfloor y\rfloor)^{\top}\text{\boldmath$\vartheta$}_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}(\lfloor y\rfloor)$ |  |  |  |  |
|  | $\text{\boldmath$a$}_{\text{Bs},M}(\lfloor y\rfloor)^{\top}\text{\boldmath$\vartheta$}_{1}$ |  | $\log(O_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$ | $\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$ | $\log(\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$
$\mathbb{R}^{+}$ | Survival Analysis |
| $\log(y)+\vartheta_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$ |  |  |  |  | Exponential AFT
|  | $\beta$ |  | log-OR |  | log-HR
| $(1,\log(y))^{\top}\text{\boldmath$\vartheta$}_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$ |  | log-normal AFT | log-logistic AFT |  | Weibull AFT
|  | $(1,\log(y))^{\top}\text{\boldmath$\vartheta$}_{1}$ |  | $\log(O_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$ | $\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$ | $\log(\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$
|  | $\beta$ |  | log-OR | AH | log-HR
$\mathbb{R}$ | Continuous Regression and Survival Analysis |
| $y\vartheta_{1}-\vartheta_{2}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$ |  | NLRM |  |  |
|  | $\vartheta_{1}^{-2}$ | variance |  |  |
|  | $(\vartheta_{2}+\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$})/\vartheta_{1}$ | mean |  |  |
| $\text{\boldmath$a$}_{\text{Bs},M}(y)^{\top}\text{\boldmath$\vartheta$}_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$ |  |  |  | Aalen AH | Cox PH
|  | $\text{\boldmath$a$}_{\text{Bs},M}(y)^{\top}\text{\boldmath$\vartheta$}_{1}$ |  | $\log(O_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$ | $\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$ | $\log(\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$
|  | $\beta$ |  | log-OR | AH | log-HR
| $\text{\boldmath$a$}_{\text{Bs},M}(y)^{\top}\text{\boldmath$\vartheta$}_{1}-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}(y)$ |  | Distribution Regression | Time-varying Cox
|  | $\text{\boldmath$a$}_{\text{Bs},M}(y)^{\top}\text{\boldmath$\vartheta$}_{1}$ |  | $\log(O_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0}))$ | $\Lambda_{Y}(y\mid\text{\boldmath$X$}=\mathbf{0})$ | $\log(\Lambda_{Y}(y\mid\text{\boldmath$x$}=\mathbf{0}))$
Table 1: Non-exhaustive Overview of Conditional Transformation Models.
Abbreviations: Proportional hazards (PH), proportional odds (PO), additive
hazards (AH), odds ratio (OR), hazard ratio (HR), complementary log (clog),
complementary log-log (cloglog), normal linear regression model (NLRM), binary
generalised linear model (BGLM), accelerated failure time (AFT).

##  5 Empirical Evaluation

We will illustrate the range of possible applications of likelihood-based
conditional transformation models. In Subsection 5.2, we will present a small
simulation experiment highlighting the possible advantage of indirectly
modelling conditional distributions with transformation functions.

###  5.1 Illustrations

##### Density Estimation: Old Faithful Geyser

The duration of eruptions and the waiting time between eruptions of the Old
Faithful geyser in the Yellowstone National Park became a standard benchmark
for non-parametric density estimation (the original data were given by
Azzalini and Bowman, 1990). The nine parameters of the transformation model
$(\Phi,\text{\boldmath$a$}_{\text{Bs},8}(\text{waiting}),\text{\boldmath$\vartheta$})$
were fitted by maximisation of the approximate log-likelihood (treating the
waiting times as “exact” observations) under the eight linear constraints
$\text{\boldmath$D$}_{9}\text{\boldmath$\vartheta$}>0$. The model depicted in
Figure 1A reproduces the classic bimodal unconditional density of waiting time
along with a kernel density estimate. It is important to note that the
transformation model was fitted likelihood-based, whereas the kernel density
estimate relied on a cross-validated bandwidth. An unconditional density
estimate for the duration of the eruptions needs to deal with censoring
because exact duration times are only available for the daytime measurements.
At night, the observations were either left-censored (“short” eruption),
interval-censored (“medium” eruption) or right-censored (“long” eruption) as
explained by Azzalini and Bowman (1990). This censoring was widely ignored in
analyses of the Old Faithful data because most non-parametric kernel
techniques cannot deal with censoring (see for example Hyndman and Yao, 2002).
We applied the transformation model
$(\Phi,\text{\boldmath$a$}_{\text{Bs},8}(\text{duration}),\text{\boldmath$\vartheta$})$
based on the exact log-likelihood function under eight linear constraints and
obtained the unconditional density depicted in Figure 1B. In Hothorn (2017b),
results for $M=40$ are computed, which led to almost identical estimates of
the distribution function.

Figure 1: Old Faithful Geyser. Estimated density for waiting times (A) between
and duration (B) of eruptions by the most likely transformation model (MLT)
and kernel smoothing (function npudens() from package np in R). Note that the
kernel estimator was based on the imputed duration times $2,3$ and $4$ for
short, medium and long eruptions at night (as are the rugs in B).

In addition, we modelled the conditional distribution of the censored duration
times given the waiting times using the transformation model
$(\Phi,(\text{\boldmath$a$}_{\text{Bs},7}(\text{duration})^{\top}\otimes\,\text{\boldmath$b$}_{\text{Bs},3}(\text{waiting})^{\top})^{\top},\text{\boldmath$\vartheta$})$.
This conditional transformation model allows a smooth conditional distribution
function of duration smoothly varying with waiting time and was fitted by
maximisation of the exact log-likelihood under $7\times 4$ linear constraints.
The corresponding conditional density in Figure 2 shows that the marginal
bimodality can be explained by relatively long durations for short waiting
times and two clusters of short and long durations after waiting times longer
than $70$ minutes. Note that in this model, the choice of $F_{Z}$ does not
influence the estimated conditional distribution or density function if the
parameterisation of $h$ is flexible enough to compensate for this change.

Figure 2: Old Faithful Geyser. Conditional distribution of duration given
waiting time estimated with a most likely transformation model (MLT) and a
conditional kernel density estimator (function npcdens() from package np).
Exact observations are given as dots, censored observations as lines.

##### Quantile Regression: Head Circumference

The Fourth Dutch Growth Study (Fredriks _et al._ , 2000) is a cross-sectional
study on growth and development of the Dutch population younger than $22$
years. Stasinopoulos and Rigby (2007) fitted a growth curve to head
circumferences (HC) of $7040$ boys using a GAMLSS model with a Box-Cox $t$
distribution describing the first four moments of head circumference
conditionally on age. The model showed evidence of kurtosis, especially for
older boys. We fitted the same growth curves by the conditional transformation
model
$(\Phi,(\text{\boldmath$a$}_{\text{Bs},3}(\text{HC})^{\top}\otimes\text{\boldmath$b$}_{\text{Bs},3}(\text{age}^{1/3})^{\top})^{\top},\text{\boldmath$\vartheta$})$
by maximisation of the approximate log-likelihood under $3\times 4$ linear
constraints
$(\text{\boldmath$D$}_{4}\otimes\text{\boldmath$I$}_{4})\text{\boldmath$\vartheta$}>0$.
Figure 3 shows the data overlaid with quantile curves obtained via inversion
of the estimated conditional distributions. The figure very closely reproduces
the growth curves presented in Figure 16 of Stasinopoulos and Rigby (2007) and
also indicates a certain asymmetry towards older boys.

Figure 3: Head Circumference Growth. Observed head circumference and age of
$7040$ boys with estimated quantile curves for
$p=0.04,0.02,0.1,0.25,0.5,0.75,0.9,0.98,0.996$.

##### Survival Analysis: German Breast Cancer Study Group-2 Trial

This prospective, controlled clinical trial on the treatment of node-positive
breast cancer patients was conducted by the German Breast Cancer Study Group
(GBSG-2, Schumacher _et al._ , 1994). Patients not older than $65$ years with
positive regional lymph nodes but no distant metastases were included in the
study. Out of $686$ women, $246$ received hormonal therapy, whereas the
control group of $440$ women did not. Additional variables include age,
menopausal status, tumour size, tumour grade, number of positive lymph nodes,
progesterone receptor and oestrogen receptor. The right-censored recurrence-
free survival time is the response variable of interest.

The Cox model
$(F_{\operatorname{MEV}},(\text{\boldmath$a$}_{\text{Bs},10}^{\top},\mathds{1}(\text{hormonal
therapy}))^{\top},\text{\boldmath$\vartheta$})$ implements the transformation
function
$h(y\mid\text{treatment})=\text{\boldmath$a$}_{\text{Bs},10}(y)^{\top}\text{\boldmath$\vartheta$}_{1}+\mathds{1}(\text{hormonal
therapy})\beta$, where
$\text{\boldmath$a$}_{\text{Bs},10}^{\top}\text{\boldmath$\vartheta$}_{1}$ is
the log-cumulative baseline hazard function parameterised by a Bernstein
polynomial and $\beta\in\mathbb{R}$ is the log-hazard ratio of hormonal
therapy. This is the classical Cox model with one treatment parameter $\beta$
but with fully parameterised baseline transformation function, which was
fitted by the exact log-likelihood under ten linear constraints. The model
assumes proportional hazards, an assumption whose appropriateness we wanted to
assess using the non-proportional hazards model
$(F_{\operatorname{MEV}},(\text{\boldmath$a$}_{\text{Bs},10}^{\top}\otimes(1,\mathds{1}(\text{hormonal
therapy})))^{\top},\text{\boldmath$\vartheta$})$ with the transformation
function

$\displaystyle
h(y\mid\text{treatment})=\text{\boldmath$a$}_{\text{Bs},10}(y)^{\top}\text{\boldmath$\vartheta$}_{1}+\mathds{1}(\text{hormonal
therapy})\text{\boldmath$a$}_{\text{Bs},10}(y)^{\top}\text{\boldmath$\vartheta$}_{2}.$

The function
$\text{\boldmath$a$}_{\text{Bs},10}(y)^{\top}\text{\boldmath$\vartheta$}_{2}$
is the time-varying difference of the log-hazard functions of women without
and with hormonal therapy and can be interpreted as the deviation from a
constant log-hazard ratio treatment effect of hormonal therapy. Under the null
hypothesis of no treatment effect, we would expect
$\text{\boldmath$\vartheta$}_{2}\equiv\mathbb{0}$. This monotonic deviation
function adds ten linear constraints
$\text{\boldmath$D$}_{11}\text{\boldmath$\vartheta$}_{1}+\text{\boldmath$D$}_{11}\text{\boldmath$\vartheta$}_{2}>0$,
which also ensure monotonicity of the transformation function for treated
patients. We first compared the fitted survivor functions obtained from the
model including a time-varying treatment effect with the Kaplan-Meier
estimators in both treatment groups. Figure 4A shows a nicely smoothed version
of the survivor functions obtained from this transformation model. Figure 4B
shows the time-varying treatment effect
$\text{\boldmath$a$}_{\text{Bs},10}(y)^{\top}\hat{\text{\boldmath$\vartheta$}}_{2}$,
together with a $95\%$ confidence band computed from the joint normal
distribution of $\hat{\text{\boldmath$\vartheta$}}_{2}$ for a grid over time
as described by Hothorn _et al._ (2008); the method is much simpler than other
methods for inference on time-varying effects (e.g. Sun _et al._ , 2009). The
$95\%$ confidence interval around the log-hazard ratio $\hat{\beta}$ is also
plotted, and as the latter is fully covered by the confidence band for the
time-varying treatment effect, there is no reason to question the treatment
effect computed under the proportional hazards assumption. An alternative
method for this type of analysis has been recently suggested by Yang and
Prentice (2015).

Figure 4: GBSG-2. Estimated survivor functions by the most likely
transformation model (MLT) and the Kaplan-Meier (KM) estimator in the two
treatment groups (A). Verification of proportional hazards (B): The log-hazard
ratio $\hat{\beta}$ (dashed line) with $95\%$ confidence interval (dark grey)
is fully covered by a $95\%$ confidence band for the time-varying treatment
effect (the time-varying log-hazard ratio is in light grey, the estimate is
the solid line) computed from a non-proportional hazards model.

In the second step, we allowed an age-varying treatment effect to be included
in the model
$(F_{\operatorname{MEV}},(\text{\boldmath$a$}_{\text{Bs},10}(y)^{\top}\otimes(\mathds{1}(\text{hormonal
therapy}),1-\mathds{1}(\text{hormonal
therapy}))\otimes\text{\boldmath$b$}_{\text{Bs},3}(\text{age})^{\top})^{\top},\text{\boldmath$\vartheta$})$.
For both treatment groups, we estimated a conditional transformation function
of survival time $y$ given age parameterised as the tensor basis of two
Bernstein bases. Each of the two basis functions comes with $10\times 3$
linear constraints; therefore, the model was fitted under $60$ linear
constraints. Figure 5 allows an assessment of the prognostic and predictive
properties of age. As the survivor functions were clearly larger for all
patients treated with hormones, the positive treatment effect applied to all
patients. However, the size of the treatment effect varied greatly. The effect
was most pronounced for women younger than $30$ and levelled off a little for
older patients. In general, the survival times were longest for women between
$40$ and $60$ years old. Younger women suffered the highest risk; for women
older than $60$ years, the risk started to increase again. This effect was
shifted towards younger women when hormonal treatment was applied.

Figure 5: GBSG-2. Prognostic and predictive effect of age. The contours depict
the conditional survivor functions given treatment and age of the patient.

##### Count Regression: Tree Pipit Counts

Müller and Hothorn (2004) reported data on the number of tree pipits Anthus
trivialis, a small passerine bird, counted on $86$ forest plots in a light
gradient ranging from open and sunny stands (small cover storey) to dense and
dark stands (large cover storey). We modelled the conditional distribution of
the number of tree pipits on one plot given the cover storey on this plot by
the transformation model
$(\Phi,(\text{\boldmath$a$}^{\top}\otimes\text{\boldmath$b$}_{\text{Bs},4}(\text{cover
storey})^{\top})^{\top},\text{\boldmath$\vartheta$})$, where
$\text{\boldmath$a$}(y)=\text{\boldmath$e$}_{5}(y+1),y=0,\dots,4$; the model
was fitted under $4\times 5$ linear constraints. In this model of count data,
the conditional distribution depends on both the number of counted birds and
the cover storey and the effect of cover storey may change with different
numbers of birds observed. Figure 6A depicts the observations, and Figure 6B
shows the conditional distribution function evaluated for $0,\dots,5$ observed
birds. The conditional distribution function obtained from a generalised
additive Poisson (GAM) model with smooth mean effect of cover storey is given
in Figure 6C. Despite some overfitting, this model is more restrictive than
our transformation model because one mean function determines the whole
distribution (the local minima of the conditional distributions as a function
of cover storey were constant in Figure 6C, whereas they were shifted towards
higher values of cover storey in Figure 6B).

Figure 6: Tree Pipit Counts. Observations (A, the size of the points is
proportional to the number of observations) and estimated conditional
distribution of number of tree pipits given cover storey by the most likely
transformation model (MLT, B) and a generalised additive Poisson model
(function gam() in package mgcv, GAM, C).

###  5.2 Simulation Experiment

The transformation family includes linear as well as very flexible models, and
we therefore illustrate the potential gain of modelling a transformation
function $h$ by comparing a very simple transformation model to a fully
parametric approach and to a non-parametric approach using a data-generating
process introduced by Hothorn _et al._ (2014).

In the transformation model
$(\Phi,((1,y)\otimes(1,\text{\boldmath$x$}^{\top}))^{\top},\text{\boldmath$\vartheta$})$,
two explanatory variables $\text{\boldmath$x$}=(x_{1},x_{2})^{\top}$ influence
both the conditional mean and the conditional variance of a normal response
$Y$. Although the transformation function is linear in $y$ with three linear
constraints, the mean and variance of $Y$ given $x$ depend on $x$ in a non-
linear way. The choices
$x_{1}\sim\operatorname{U}[0,1],x_{2}\sim\operatorname{U}[-2,2]$ with
$\text{\boldmath$\vartheta$}=(0,0,-1,.5,1,0)$ lead to the heteroscedastic
varying coefficient model

$\displaystyle
Y=\frac{1}{x_{1}+0.5}x_{2}+\frac{1}{x_{1}+0.5}\varepsilon,\quad\varepsilon\sim\operatorname{N}(0,1),$
(9)

where the variance of $Y$ ranges between $0.44$ and $4$ depending on $x_{1}$.
This model can be fitted in the GAMLSS framework under the assumptions that
the mean of the normal response depends on a smoothly varying regression
coefficient $(x_{1}+0.5)^{-1}$ for $x_{2}$ and that the variance is a smooth
function of $x_{1}$. This model is therefore fully parametric. As a non-
parametric counterpart, we used a kernel estimator for estimating the
conditional distribution function of $Y$ as a function of the two explanatory
variables.

From the transformation model, the GAMLSS and kernel estimators, we obtained
estimates of $F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)$ over a grid
on $y,x_{1},x_{2}$ and computed the mean absolute deviation (MAD) of the true
and estimated probabilities

$\displaystyle\text{MAD}(x_{1},x_{2})=\frac{1}{n}\sum_{y}|F_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}}(y)-\hat{F}_{Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},N}(y)|$

for each pair of $x_{1}$ and $x_{2}$. Then, the minimum, the median and the
maximum of the MAD values for all $x_{1}$ and $x_{2}$ were computed as summary
statistics. The most likely transformation approach and its two competitors
were estimated and evaluated for $100$ random samples of size $N=200$ drawn
from model (9). Cross-validation was used to determine the bandwidths for the
kernel-based estimators (function npcdist() in package np; for details, see
Hayfield and Racine, 2008). We fitted the GAMLSS models by boosting; the
number of boosting iterations was determined via sample splitting (Mayr _et
al._ , 2012). To investigate the stability of the three procedures under non-
informative explanatory variables, we added to the data $p=1,\dots,5$
uniformly distributed variables without association to the response and
included them as potential explanatory variables in the three models. The case
$p=0$ corresponds to model (9).

Figure 7: Empirical Evaluation. Minimum, median and maximum of the mean
absolute deviation (MAD) between true and estimated probabilities for most
likely transformation models (MLT), non-parametric kernel distribution
function estimation (Kernel) and generalised additive models for location,
scale and shape (GAMLSS) for $100$ random samples. Values on the ordinate can
be interpreted as absolute differences of probabilities. The grey horizontal
lines correspond to the median of MLT.

Figure 7 shows the empirical distributions of the minimum, median and maximum
MAD for the three competitors. Except for the minimum MAD in the absence of
any irrelevant explanatory variables ($p=0$), the conditional distributions
fitted by the transformation models were closer to the true conditional
distribution function by means of the MAD. This result was obtained because
the transformation model only had to estimate a simple transformation
function, whereas the other two procedures had a difficult time approximating
this simple transformation model on another scale. However, the comparison
illustrates the potential improvement one can achieve when fitting simple
models for the transformation function instead of more complex models for the
mean (GAMLSS) or distribution function (Kernel). The kernel estimator led to
the largest median MAD values but seemed more robust than GAMLSS with respect
to the maximum MAD. These results were remarkably robust in the presence of up
to five non-informative explanatory variables, although of course the MAD
increased with the number of non-informative variables $p$.

##  6 Discussion

The contribution of a likelihood approach for the general class of conditional
transformation models is interesting both from a theoretical and a practical
perspective. With the range of simple to very complex transformation functions
introduced in Section 4 and illustrated in Section 5, it becomes possible to
understand classical parametric, semi-parametric and non-parametric models as
special cases of the same model class. Thus, analytic comparisons between
models of different complexity become possible. The transformation family
$\mathbb{P}_{Y,\Theta}$, the corresponding likelihood function and the most
likely transformation estimator are easy to understand. This makes the
approach appealing also from a teaching perspective. Connections between
standard parametric models (for example, the normal linear model) and
potentially complex models for survival or ordinal data can be outlined in
very simple notation, placing emphasis on the modelling of (conditional)
distributions instead of just modelling (conditional) means. Computationally,
the log-likelihood $\log\circ\mathcal{L}$ is linear in the number of
observations $N$ and, for contributions of “exact continuous” responses, only
requires the evaluation of the derivative $h^{\prime}$ of the transformation
function $h$ instead of integrals thereof. Standard optimisers for linearly
constrained problems can be applied, and AIC- or BIC-based model selection is
possible. Transformation models directly suggest a way to specify conceptual
$(F_{Z},\text{\boldmath$c$},\text{\boldmath$\vartheta$})$ and fitted
$(F_{Z},\text{\boldmath$c$},\hat{\text{\boldmath$\vartheta$}})$ models in
computer languages because “only” the conditional transformation function
$h=\text{\boldmath$c$}^{\top}\text{\boldmath$\vartheta$}$ has to be specified
in addition to a “simple” distribution function $F_{Z}$ (Hothorn, 2017a, b).

The results presented in Section 5 are based on only roughly $1000$ lines of R
code. We tested this reference implementation of most likely transformations
against all linear transformation models available in R to date and obtained
equivalent regression coefficients $\hat{\text{\boldmath$\beta$}}$ and their
corresponding covariance matrices. Thus, the framework helps to reduce the
code base and thereby helps to considerably reduce the number of possible
errors in implementations of linear transformation models, while at the same
time it allows novel and more complex models to be fitted (Hothorn, 2017a, b).

Based on the general understanding of transformation models outlined in this
paper, it will be interesting to study these models outside the strict
likelihood world. A mixed transformation model for cluster data (Cai _et al._
, 2002; Huber-Carol and Vonta, 2004; Zeng _et al._ , 2005; Choi and Huang,
2012) is often based on the transformation function
$h(y\mid\text{\boldmath$x$},i)=h_{Y}(y)+\delta_{i}+h_{\text{\boldmath$x$}}(\text{\boldmath$x$})$
with random intercept (or “frailty” term) $\delta_{i}$ for the $i$th
observational unit. Conceptually, a more complex deviation from the global
model could by formulated as
$h(y\mid\text{\boldmath$x$},i)=h_{Y}(y)+h_{Y}(y,i)+h_{\text{\boldmath$x$}}(\text{\boldmath$x$})$,
i.e. each observational unit is assigned its own “baseline” transformation
$h_{Y}(y)+h_{Y}(y,i)$, where the second term is an integral zero deviation
from $h_{Y}$. For longitudinal data with possibly time-varying explanatory
variables, the model
$h(y\mid\text{\boldmath$x$}(t),t)=h_{Y}(y,t)+\text{\boldmath$x$}(t)\text{\boldmath$\beta$}(t)$
(Ding _et al._ , 2012; Wu and Tian, 2013) can also be understood as a mixed
version of a conditional transformation model. The penalised log-likelihood
$\log(\mathcal{L}(h\mid y))-\text{pen}(\text{\boldmath$\beta$})$ for the
linear transformation model
$h(y\mid\text{\boldmath$x$})=h_{Y}(y)-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}$
leads to Ridge- or Lasso-type regularised models, depending on the form of the
penalty term. Priors for all model parameters $\vartheta$ allow a fully
Bayesian treatment of transformation models. Instead of the relatively strict
model assumption
$F_{Y}=F_{Z}\circ\text{\boldmath$a$}^{\top}\text{\boldmath$\vartheta$}$, one
could of course allow arbitrary unknown transformation functions
$h\in\mathcal{H}$ in the spirit of non-parametric regression and study the
quality of the approximation
$\text{\boldmath$a$}^{\top}\hat{\text{\boldmath$\vartheta$}}$ for $h$ for the
fixed and random design case.

It is possible to relax the assumption that $F_{Z}$ is known. The simultaneous
estimation of $F_{Z}$ in the model $\mathbb{P}(Y\leq
y\mid\text{\boldmath$X$}=\text{\boldmath$x$})=F_{Z}(h_{Y}(y)-\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$})$
was studied by Horowitz (1996) and later extended by Linton _et al._ (2008) to
non-linear functions $h_{\text{\boldmath$x$}}$ with parametric baseline
transformation $h_{Y}$ and kernel estimates for $F_{Z}$ and
$h_{\text{\boldmath$x$}}$. For AFT models, Zhang and Davidian (2008) applied
smooth approximations for the density $f_{Z}$ in an exact censored likelihood
estimation procedure. In a similar setup, Huang (2014) proposed a method to
jointly estimate the mean function and the error distribution in a generalised
linear model. The estimation of $F_{Z}$ is noteworthy in additive models of
the form $h_{Y}+h_{\text{\boldmath$x$}}$ because these models assume
additivity of the contributions of $y$ and $x$ on the scale of
$F_{Z}^{-1}(\mathbb{P}(Y\leq y\mid\text{\boldmath$X$}=\text{\boldmath$x$}))$.
If this model assumption seems questionable, one can either allow unknown
$F_{Z}$ or move to a transformation model featuring a more complex
transformation function. From this point of view, the distribution function
$F_{Z}$ in flexible transformation models is only a computational device
mapping the unbounded transformation function $h$ into the unit interval
strictly monotonically, making the evaluation of the likelihood easy. Then,
$F_{Z}$ has no further meaning or interpretation as error distribution. A
compromise could be the family of distributions
$F_{Z}(z\mid\rho)=1-(1+\rho\exp(z))^{-1/\rho}$ for $\rho>0$ (suggested by
McLain and Ghosh, 2013) with simultaneous maximum likelihood estimation of
$\vartheta$ and $\rho$ for additive transformation functions
$h=h_{Y}+h_{\text{\boldmath$x$}}$, as these models are flexible and still
relatively easy to interpret.

Another comment concerns the connection between transformation models and
quantile regression. For some probability $p\in[0,1]$, the conditional
quantile function is

$\displaystyle
Q_{Y}(p\mid\text{\boldmath$X$}=\text{\boldmath$x$}):=F_{Y}^{-1}(p\mid\text{\boldmath$X$}=\text{\boldmath$x$})=\inf\\{y\in\Xi\mid
F_{Z}(h(y\mid\text{\boldmath$x$}))\geq p\\},$

and for absolutely continuous $Y$, we then get
$Q_{Y}(p\mid\text{\boldmath$X$}=\text{\boldmath$x$})=h^{-1}(F_{Z}^{-1}(p)\mid\text{\boldmath$X$}=\text{\boldmath$x$})$.
A linear quantile regression model (Koenker, 2005) assumes
$Q_{Y}(p\mid\text{\boldmath$X$}=\text{\boldmath$x$})=\alpha(p)+\tilde{\text{\boldmath$x$}}^{\top}\text{\boldmath$\beta$}(p)$,
which cannot be written as a linear transformation model. Thus, non-linear
transformation functions are necessary in order to achieve the same
flexibility as a linear quantile regression model; however, the full
conditional distribution function can be estimated in one step, which avoids
computational problems such as quantile crossing. In the most complex case of
a transformation model that is invariant with respect to the choice of
$F_{Z}$, there is a correspondingly complex quantile regression model such
that both models are equivalent, as was also noted by Chernozhukov _et al._
(2013). In a certain sense, we can understand transformation models as
“inverse quantile regression”. A more detailed analysis of the connection
between transformation models (called distribution regression there) and
quantile regression can be found in Leorato and Peracchi (2015).

Figure 8: GBSG-2. Comparison of exact and partial likelihood for order
$M=1,\dots,30,35,40,45,50$ of the Bernstein polynomial approximating the log-
cumulative baseline hazard function $h_{Y}$. The estimated log-cumulative
baseline hazard functions for varying $M$ are shown in grey and the Nelson-
Aalen-Breslow estimator is shown in red (A). The right panel (B) shows the
trajectories of the regression coefficients $\hat{\text{\boldmath$\beta$}}$
obtained for varying $M$, which are represented as dots. The horizontal lines
represent the partial likelihood estimates.

In light of the empirical results discussed in this paper and the theoretical
work of McLain and Ghosh (2013) on a Cox model with log-cumulative baseline
hazard function parameterised in terms of a Bernstein polynomial with
increasing order $M$, one might ask where the boundaries between parametric,
semi-parametric and non-parametric statistics lie. The question how the order
$M$ affects results practically has been repeatedly raised; therefore, we will
close our discussion by looking at a Cox model with increasing $M$ for the
GBSG-2 data. All eight baseline variables were included in the linear
predictor, and we fitted the model with orders $M=1,\dots,30,35,40,45,50$ of
the Bernstein polynomial parameterising the log-cumulative baseline hazard
function. In Figure 8A, the log-cumulative baseline hazard functions start
with a linear function ($M=1$) and quickly approach a function that is
essentially a smoothed version of the Nelson-Aalen-Breslow estimator plotted
in red. In Figure 8B, the trajectories of the estimated regression
coefficients become very similar to the partial likelihood estimates as $M$
increased. For $M\geq 10$, for instance, the results of the “semi-parametric”
and the “fully parametric” Cox models are practically equivalent. An extensive
collection of such head-to-head comparisons of most likely transformations
with their classical counterparts can be found in Hothorn (2017b). Our work
for this paper and practical experience with its reference software
implementation convinced us that rethinking classical models in terms of fully
parametric transformations is intellectually and practically a fruitful
exercise.

## Acknowledgements

Torsten Hothorn received financial support by Deutsche Forschungsgemeinschaft
(DFG) under grant number HO 3242/4-1. We thank Karen A. Brune for improving
the language.

## References

  * Alzaatreh _et al._ (2013) Alzaatreh A, Lee C, Famoye F (2013).  “A New Method for Generating Families of Continuous Distributions.”  _Metron_ , 71(1), 63–79.  10.1007/s40300-013-0007-y.
  * Aranda-Ordaz (1983) Aranda-Ordaz FJ (1983).  “An Extension of the proportional Hazards Model for Grouped Data.”  _Biometrics_ , 39(1), 109–117.
  * Azzalini and Bowman (1990) Azzalini A, Bowman AW (1990).  “A Look at Some Data on the Old Faithful Geyser.”  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 39(3), 357–365.  10.2307/2347385.
  * Babu _et al._ (2002) Babu GJ, Canty AJ, Chaubey YP (2002).  “Application of Bernstein Polynomials for Smooth Estimation of a Distribution and Density Function.”  _Journal of Statistical Planning and Inference_ , 105(2), 377–392.  10.1016/S0378-3758(01)00265-8.
  * Bender _et al._ (2005) Bender R, Augustin T, Blettner M (2005).  “Generating Survival Times to Simulate Cox Proportional Hazards Models.”  _Statistics in Medicine_ , 24(11), 1713–1723.  10.1002/sim.2059.
  * Bickel _et al._ (1993) Bickel PJ, Klaassen CAJ, Ritov Y, Wellner JA (1993).  _Efficient and Adaptive Estimation for Semiparametric Models_.  The Johns Hopkins University Press, Baltimore, U.S.A. and London, U.K.
  * Box and Cox (1964) Box GEP, Cox DR (1964).  “An Analysis of Transformations.”  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 26(2), 211–252.
  * Cai and Betensky (2003) Cai T, Betensky RA (2003).  “Hazard Regression for Interval-Censored Data with Penalized Spline.”  _Biometrics_ , 59(3), 570–579.  10.1111/1541-0420.00067.
  * Cai _et al._ (2002) Cai T, Cheng SC, Wei LJ (2002).  “Semiparametric Mixed-Effects Models for Clustered Failure Time Data.”  _Journal of the American Statistical Association_ , 97(458), 514–522.  10.1198/016214502760047041.
  * Chang _et al._ (2005) Chang I, Hsiung CA, Wu YJ, Yany CC (2005).  “Bayesian Survival Analysis Using Bernstein Polynomials.”  _Scandinavian Journal of Statistics_ , 32(3), 447–466.  10.1111/j.1467-9469.2005.00451.x.
  * Chen _et al._ (2002) Chen K, Jin Z, Ying Z (2002).  “Semiparametric Analysis of Transformation Models with Censored Data.”  _Biometrika_ , 89(3), 659–668.  10.1093/biomet/89.3.659.
  * Cheng _et al._ (1995) Cheng SC, Wei LJ, Ying Z (1995).  “Analysis of Transformation Models with Censored Data.”  _Biometrika_ , 82(4), 835–845.  10.1093/biomet/82.4.835.
  * Cheng _et al._ (1997) Cheng SC, Wei LJ, Ying Z (1997).  “Predicting Survival Probabilities with Semiparametric Transformation Models.”  _Journal of the American Statistical Association_ , 92(437), 227–235.  10.1080/01621459.1997.10473620.
  * Chernozhukov _et al._ (2013) Chernozhukov V, Fernández-Val I, Melly B (2013).  “Inference on Counterfactual Distributions.”  _Econometrica_ , 81(6), 2205–2268.  10.3982/ECTA10582.
  * Choi and Huang (2012) Choi S, Huang X (2012).  “A General Class of Semiparametric Transformation Frailty Models for Nonproportional Hazards Survival Data.”  _Biometrics_ , 68(4), 1126–1135.  10.1111/j.1541-0420.2012.01784.x.
  * Cox (1975) Cox DR (1975).  “Partial Likelihood.”  _Biometrika_ , 62(2), 269–276.
  * Crowther and Lambert (2014) Crowther MJ, Lambert PC (2014).  “A General Framework for Parametric Survival Analysis.”  _Statistics in Medicine_ , 33(30), 5280–5297.  10.1002/sim.6300.
  * Curtis and Ghosh (2011) Curtis SM, Ghosh SK (2011).  “A Variable Selection Approach to Monotonic Regression with Bernstein Polynomials.”  _Journal of Applied Statistics_ , 38(5), 961–976.  10.1080/02664761003692423.
  * Ding _et al._ (2012) Ding AA, Tian S, Yu Y, Guo H (2012).  “A Class of Discrete Transformation Survival Models With Application to Default Probability Prediction.”  _Journal of the American Statistical Association_ , 107(499), 990–1003.  10.1080/01621459.2012.682806.
  * Doksum and Gasko (1990) Doksum KA, Gasko M (1990).  “On a Correspondence Between Models in Binary Regression Analysis and in Survival Analysis.”  _International Statistical Review_ , 58(3), 243–252.
  * Farouki (2012) Farouki RT (2012).  “The Bernstein Polynomial Basis: A Centennial Retrospective.”  _Computer Aided Geometric Design_ , 29(6), 379–419.  10.1016/j.cagd.2012.03.001.
  * Fisher (1934) Fisher RA (1934).  “Two New Properties of Mathematical Likelihood.”  _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ , 53, 285–306.
  * Foresi and Peracchi (1995) Foresi S, Peracchi F (1995).  “The Conditional Distribution of Excess Returns: An Empirical Analysis.”  _Journal of the American Statistical Association_ , 90(430), 451–466.  10.1080/01621459.1995.10476537.
  * Fraser (1968) Fraser DAS (1968).  _The Structure of Inference_.  John Wiley & Sons, New York, U.S.A.
  * Fredriks _et al._ (2000) Fredriks AM, van Buuren S, Burgmeijer RJF, Meulmeester JF, Beuker RJ, Brugman E, Roede MJ, Verloove-Vanhorick SP, Wit J (2000).  “Continuing Positive Secular Growth Change in The Netherlands 1955–1997.”  _Pediatric Research_ , 47(3), 316–323.
  * Gneiting and Katzfuss (2014) Gneiting T, Katzfuss M (2014).  “Probabilistic Forecasting.”  _Annual Review of Statistics and Its Application_ , 1(1), 125–151.  10.1146/annurev-statistics-062713-085831.
  * Hayfield and Racine (2008) Hayfield T, Racine JS (2008).  “Nonparametric Econometrics: The np Package.”  _Journal of Statistical Software_ , 27(5), 1–32.  URL http://www.jstatsoft.org/v27/i05.
  * Hofner _et al._ (2016) Hofner B, Mayr A, Fenske N, Schmid M (2016).  _gamboostLSS: Boosting Methods for ’GAMLSS’_.  R package version 1.2-2, URL https://CRAN.R-project.org/package=gamboostLSS.
  * Horowitz (1996) Horowitz JL (1996).  “Semiparametric Estimation of a Regression Model with an Unknown Transformation of the Dependent Variable.”  _Econometrica_ , 64(1), 103–137.
  * Hothorn (2017a) Hothorn T (2017a).  _mlt: Most Likely Transformations_.  R package version 0.1-3, URL https://CRAN.R-project.org/package=mlt.
  * Hothorn (2017b) Hothorn T (2017b).  _Most Likely Transformations: The mlt Package_.  R package vignette version 0.1-5, URL https://CRAN.R-project.org/package=mlt.docreg.
  * Hothorn _et al._ (2008) Hothorn T, Bretz F, Westfall P (2008).  “Simultaneous Inference in General Parametric Models.”  _Biometrical Journal_ , 50(3), 346–363.  10.1002/bimj.200810425.
  * Hothorn _et al._ (2014) Hothorn T, Kneib T, Bühlmann P (2014).  “Conditional Transformation Models.”  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(1), 3–27.  10.1111/rssb.12017.
  * Huang (2014) Huang A (2014).  “Joint Estimation of the Mean and Error Distribution in Generalized Linear Models.”  _Journal of the American Statistical Association_ , 109(505), 186–196.  10.1080/01621459.2013.824892.
  * Huber-Carol and Vonta (2004) Huber-Carol C, Vonta I (2004).  “Frailty Models for Arbitrarily Censored and Truncated Data.”  _Lifetime Data Analysis_ , 10(4), 369–388.  10.1007/s10985-004-4773-y.
  * Hyndman and Yao (2002) Hyndman RJ, Yao Q (2002).  “Nonparametric Estimation and Symmetry Tests for Conditional Density Functions.”  _Journal of Nonparametric Statistics_ , 14(3), 259–278.  10.1080/10485250212374.
  * Joly _et al._ (1998) Joly P, Commenges D, Letenneur L (1998).  “A Penalized Likelihood Approach for Arbitrarily Censored and Truncated Data: Application to Age-Specific Incidence of Dementia.”  _Biometrics_ , 54(1), 185–194.
  * Jones and Pewsey (2009) Jones MC, Pewsey A (2009).  “Sinh-arcsinh Distributions.”  _Biometrika_ , 96(4), 761–780.  10.1093/biomet/asp053.
  * Klein and Moeschberger (2003) Klein JP, Moeschberger MK (2003).  _Survival Analysis. Techniques for Censored and Truncated Data._ 2nd edition. Springer, New York, U.S.A.
  * Koenker (2005) Koenker R (2005).  _Quantile Regression_.  Economic Society Monographs. Cambridge University Press, New York, U.S.A.
  * Kooperberg and Clarkson (1997) Kooperberg C, Clarkson DB (1997).  “Hazard Regression with Interval-Censored Data.”  _Biometrics_ , 53(4), 1485–1494.
  * Kooperberg _et al._ (1995) Kooperberg C, Stone CJ, Truong YK (1995).  “Hazard Regression.”  _Journal of the American Statistical Association_ , 90(429), 78–94.
  * Lehmann (1983) Lehmann EL (1983).  _Theory of Point Estimation_.  John Wiley & Sons, New York, U.S.A.
  * Leorato and Peracchi (2015) Leorato S, Peracchi F (2015).  “Comparing Distribution and Quantile Regression.”  _Technical Report 1511_ , Einaudi Institute for Economics and Finance, Rome, Italy.  URL https://ideas.repec.org/p/eie/wpaper/1511.html.
  * Lindsey (1996) Lindsey JK (1996).  _Parametric Statistical Inference_.  Clarendon Press, Oxford, UK.
  * Lindsey (1999) Lindsey JK (1999).  “Some Statistical Heresies.”  _Journal of the Royal Statistical Society: Series D (The Statistician)_ , 48(1), 1–40.
  * Linton _et al._ (2008) Linton O, Sperlich S, van Keilegom I (2008).  “Estimation of a Semiparametric Transformation Model.”  _The Annals of Statistics_ , 36(2), 686–718.  10.1214/009053607000000848.
  * Ma _et al._ (2014) Ma J, Heritier S, Lô SN (2014).  “On the Maximum Penalized Likelihood Approach for Proportional Hazard Models with Right Censored Survival Data.”  _Computational Statistics & Data Analysis_, 74, 142–156.  10.1016/j.csda.2014.01.005.
  * Mallick and Walker (2003) Mallick BK, Walker S (2003).  “A Bayesian Semiparametric Transformation Model Incorporating Frailties.”  _Journal of Statistical Planning and Inference_ , 112(1–2), 159–174.  10.1016/S0378-3758(02)00330-0.
  * Mayr _et al._ (2012) Mayr A, Fenske N, Hofner B, Kneib T, Schmid M (2012).  “GAMLSS for High-dimensional Data–A Flexible Approach Based on Boosting.”  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 61(3), 403–427.  10.1111/j.1467-9876.2011.01033.x.
  * McLain and Ghosh (2013) McLain AC, Ghosh SK (2013).  “Efficient Sieve Maximum Likelihood Estimation of Time-Transformation Models.”  _Journal of Statistical Theory and Practice_ , 7(2), 285–303.
  * Möst and Hothorn (2015) Möst L, Hothorn T (2015).  “Conditional Transformation Models for Survivor Function Estimation.”  _International Journal of Biostatistics_ , 11(1), 23–50.  10.1515/ijb-2014-0006.
  * Möst _et al._ (2016) Möst L, Schmid M, Faschingbauer F, Hothorn T (2016).  “Predicting Birth Weight with Conditionally Linear Transformation Models.”  _Statistical Methods in Medical Research_ , 25(6), 2781–2810.  10.1177/0962280214532745.
  * Müller and Hothorn (2004) Müller J, Hothorn T (2004).  “Maximally Selected Two-Sample Statistics as a new Tool for the Identification and Assessment of Habitat Factors with an Application to Breeding Bird Communities in Oak Forests.”  _European Journal of Forest Research_ , 123, 218–228.  10.1007/s10342-004-0035-5.
  * Pitman (1939) Pitman EJG (1939).  “The Estimation of the Location and Scale Parameters of a Continuous Population of Any Given Form.”  _Biometrika_ , 30(3/4), 391–421.
  * R Core Team (2017) R Core Team (2017).  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria.  URL http://www.R-project.org/.
  * Racine and Hayfield (2014) Racine JS, Hayfield T (2014).  _np: Nonparametric kernel smoothing methods for mixed data types_.  R package version 0.60-2, URL https://CRAN.R-project.org/package=np.
  * Rigby and Stasinopoulos (2004) Rigby RA, Stasinopoulos DM (2004).  “Smooth Centile Curves for Skew and Kurtotic Data Modelled Using the Box-Cox Power Exponential Distribution.”  _Statistics in Medicine_ , 23(19), 3053–3076.  10.1002/sim.1861.
  * Rigby and Stasinopoulos (2005) Rigby RA, Stasinopoulos DM (2005).  “Generalized Additive Models for Location, Scale and Shape.”  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 54(3), 507–554.  10.1111/j.1467-9876.2005.00510.x.
  * Rothe and Wied (2013) Rothe C, Wied D (2013).  “Misspecification Testing in a Class of Conditional Distributional Models.”  _Journal of the American Statistical Association_ , 108(501), 314–324.  10.1080/01621459.2012.736903.
  * Royston and Parmar (2002) Royston P, Parmar MKB (2002).  “Flexible Parametric Proportional-hazards and Proportional-odds Models for Censored Survival Data, with Application to Prognostic Modelling and Estimation of Treatment Effects.”  _Statistics in Medicine_ , 21(15), 2175–2197.  10.1002/sim.1203.
  * Sabanés Bové and Held (2013) Sabanés Bové D, Held L (2013).  “Comment on Cai and Betensky (2003), On the Poisson Approximation for Hazard Regression.”  _Biometrics_ , 69(3), 795–795.  10.1111/biom.12088.
  * Schumacher _et al._ (1994) Schumacher M, Basert G, Bojar H, Hübner K, Olschewski M, Sauerbrei W, Schmoor C, Beyerle C, Neumann RLA, Rauschecker, H F, for the German Breast Cancer Study Group (1994).  “Randomized $2\times 2$ trial evaluating hormonal treatment and the duration of chemotherapy in node-positive breast cancer patients.”  _Journal of Clinical Oncology_ , 12, 2086–2093.
  * Stasinopoulos and Rigby (2007) Stasinopoulos DM, Rigby RA (2007).  “Generalized Additive Models for Location Scale and Shape (GAMLSS) in R.”  _Journal of Statistical Software_ , 23(7), 1–46.  URL http://www.jstatsoft.org/v23/i07.
  * Sun _et al._ (2009) Sun Y, Sundaram R, Zhao Y (2009).  “Empirical Likelihood Inference for the Cox Model with Time-dependent Coefficients via Local Partial Likelihood.”  _Scandinavian Journal of Statistics_ , 36(3), 444–462.  10.1111/j.1467-9469.2008.00634.x.
  * Thas _et al._ (2012) Thas O, Neve JD, Clement L, Ottoy JP (2012).  “Probabilistic Index Models.”  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 74(4), 623–671.  10.1111/j.1467-9868.2011.01020.x.
  * Tobin (1958) Tobin J (1958).  “Estimation of Relationships for Limited Dependent Variables.”  _Econometrica_ , 26(1), 24–36.  10.2307/1907382.
  * Tutz (2012) Tutz G (2012).  _Regression for Categorical Data_.  Cambridge University Press, New York, U.S.A.
  * van de Geer (2000) van de Geer S (2000).  _Empirical Processes in M-Estimation_.  Cambridge University Press, Cambridge, UK.
  * van der Vaart (1998) van der Vaart AW (1998).  _Asymptotic Statistics_.  Cambridge University Press, Cambridge, UK.
  * Varadhan (2015) Varadhan R (2015).  _alabama: Constrained Nonlinear Optimization_.  R package version 2015.3-1, URL https://CRAN.R-project.org/package=alabama.
  * Wood (2017) Wood S (2017).  _mgcv: Mixed GAM Computation Vehicle with GCV/AIC/REML Smoothness Estimation_.  R package version 1.8-17, URL https://CRAN.R-project.org/package=mgcv.
  * Wu and Tian (2013) Wu CO, Tian X (2013).  “Nonparametric Estimation of Conditional Distributions and Rank-Tracking Probabilities With Time-Varying Transformation Models in Longitudinal Studies.”  _Journal of the American Statistical Association_ , 108(503), 971–982.
  * Yang and Prentice (2015) Yang S, Prentice RL (2015).  “Assessing Potentially Time-dependent Treatment Effect From Clinical Trials and Observational Studies for Survival Data, With Applications to the Women’s Health Initiative Combined Hormone Therapy Trial.”  _Statistics in Medicine_ , 34(11), 1801–1817.  10.1002/sim.6453.
  * Zeng and Lin (2007) Zeng D, Lin DY (2007).  “Maximum Likelihood Estimation in Semiparametric Regression Models with Censored Data.”  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 69(4), 507–564.  10.1111/j.1369-7412.2007.00606.x.
  * Zeng _et al._ (2005) Zeng D, Lin DY, Yin G (2005).  “Maximum Likelihood Estimation for the Proportional Odds Model with Random Effects.”  _Journal of the American Statistical Association_ , 100(470), 470–483.  10.1198/016214504000001420.
  * Zhang and Davidian (2008) Zhang M, Davidian M (2008).  “Smooth Semiparametric Regression Analysis for Arbitrarily Censored Time-to-Event Data.”  _Biometrics_ , 64(2), 567–576.  10.1111/j.1541-0420.2007.00928.x.

## Appendix

### Computational Details

A reference implementation of most likely transformation models is available
in the mlt package (Hothorn, 2017a). All data analyses can be reproduced in
the dynamic document Hothorn (2017b). Augmented Lagrangian Minimization