§ INTRODUCTION

Regression analysis aims to characterize the conditional distribution of a response
$\rY$ given a set of covariates $\rX$, thereby describing how changes in the
covariates propagate to the conditional distribution $\rY \given \rX = \rx$
[Fahrmeir et al., 2013]. In this paper, we present
[Kook et al., 2022], an R package for estimating a broad
class of distributional regression models for various types of responses (continuous,
survival, count, ordinal, binary) using tabular or non-tabular (image or text)
data or both. Package  is available from the Comprehensive R
Archive Network (CRAN) at <https://CRAN.R-project.org/package=deeptrafo>.
Due to a flexible tensorflow [Allaire and Tang, 2022] backend and mini-batch optimization,
not only scales well with non-tabular (imaging, text) data but also
big tabular data sets. Many well-known models fall into the class of transformation
models (TMs), such as normal linear regression (Lm), Cox proportional hazards
models (CoxPH), and proportional odds logistic regression <cit.>.
In the following, we review existing software for fitting these models.

Existing software packages
TMs for tabular data are implemented in tram [Hothorn et al., 2022] using mlt
[Hothorn, 2022] and fitted via maximum likelihood, relying on alabama
[Varadhan, 2022] and BB [Varadhan and Gilbert, 2019] for optimization. Package tram
provides an intuitive interface for fitting a multitude of distributional
regression models, ranging from shift and shift-scale [Siegfried et al., 2022]
to tensor-product models [Hothorn et al., 2014]. Several extensions of
transformation models exist. For instance, cotram for count TMs [Siegfried and Hothorn, 2020],
tramME for mixed effects TMs [Tamási and Hothorn, 2021], and tramnet as well as
tramvs for regularized TMs [Kook and Hothorn, 2021, 17]. Transformation boosting
machines [Hothorn, 2020] and transformation trees and random forests [Hothorn, 2021]
offer extensions to classical machine learning models.

Neural network-based transformation models
With the advent of (deep) neural networks and the routine collection of non-tabular
data, the idea to combine deep learning and distributional regression approaches was
adopted in several ways. For instance, Rügamer et al., 2020 parameterize distributional
regression models via neural networks, Sick et al., 2021 describe flexible deep transformation
models for continuous responses. Kook et al., 2022 focus on semi-structured regression
for ordinal responses, and Rügamer et al., 2021 extend the DCTM approach to
distributional autoregressive models for time series responses. Alternative approaches to
combining regression with neural networks including generalized additive models for location,
scale, and shape have been implemented in deepregression [Rügamer et al., 2021].
In this paper, we present , which unifies the above DCTM approaches in a single
R package. Combining distributional regression with neural network-based estimation
has many advantages, such as modularity (data analysts can easily use well-established
problem-specific neural network architectures), and easy handling of big datasets (through
mini-batch gradient descent with adaptive learning rates). Unlike models in deepregression,
s are distribution-free, they do not require the specification of a parametric family
of distributions for the response given covariates.

The rest of this paper is organized as follows. Section <ref> introduces
the statistical theory behind TMs and DCTMs. The inner workings of  are
described in Section <ref>, where several case studies on how to setup up,
fit, validate, and interpret DCTMs are presented. We present an application to binary
classification with tabular and text modalities, and an application to time series
modeling via autoregressive TMs [Rügamer et al., 2021]. The appendix contains
information on advanced usage of the package, how censored responses are handled
(Appendix <ref>) or how to warmstart or fix parameters of certain predictors
(Appendix <ref>). In Appendix <ref>, we demonstrate the package for
large tabular datasets and factors with many levels, which cannot be handled by
standard implementations of classical regression models.

§.§ Deep conditional transformation models

Transformation models [Hothorn et al., 2014, Hothorn et al., 2018] estimate the entire
conditional distribution $\rY \given \rX = \rx$ of a response $\rY$ with sample space $\calY$
given a realization $\rx \in \calX$ of covariates $\rX$, via

\begin{align}
\pYx(\ry) = \pZ(\h(\ry \given \rx)), \quad \ry \in \calY,
\end{align}

without committing to a particular parametric family of distributions in $\pYx$.
Here, the cumulative
distribution function is modelled as a composition of a continuous latent distribution $\pZ
: \RR \to [0, 1]$ and a transformation function $\h : \calY \times \calX \to \RR$, which is
monotonically non-decreasing in $\ry$.

For continuous responses $\h$ is continuous and for discrete responses $\h$ is discrete
(see Figure <ref>).

Depending on the choice of $\pZ$ and $\h$, TMs cover a wide range of well-known models
with varying complexity. For instance, choosing $\pZ = \Phi$ and $\h(\ry \given \rx) =
\sigma^{-1}(\ry - \alpha - \rx^\top\shiftparm)$, with standard deviation $\sigma > 0$,
and intercept $\alpha \in \RR$, is equivalent to a normal linear regression model.
Choosing $\pZ(\rz) = 1 - \exp(-\exp(\rz))$ with $\h(\ry \given \rx) = a + b \log\ry +
\rx^\top\shiftparm$, with intercept $a$ and slope $b > 0$, is equivalent to a Weibull
regression model [Hothorn et al., 2018]. More examples are given in
Section <ref>.

In semi-structured regression, we have access to $J$ input modalities $\rx_1, \dots, \rx_J$,
such as tabular data, images, or text, from which we construct structured (linear, sparse,
or smooth) or unstructured (neural network) predictors. We take shorthand in writing
$\rx_{1:J}$ to denote the collection of all predictors $\rx_1, \dots, \rx_J$. Note that these
inputs may be non-tabular, there may be a $j$ for which $\rx_j \in \calX_j \not\subseteq
\RR^d$. By $\calX \coloneqq \calX_1 \times \dots \times \calX_J$ we denote the entire input
space. In s, restrictions on the functional form of $\h$, the way predictors are
constructed based on the input data, lead to varying degrees of interpretability and
flexibility of the model. We begin with an example before introducing $\h$ in its most flexible
form. Consider a problem with a single tabular ($\rx_1 \in \calX_1 \subseteq \RR^p$) and a
single text modality ($\rx_2 \in \calX_2 \not\subseteq \RR^d)$. Data analysts commonly assume
additivity in the effects the separate modalities, which can be realized by modelling the effect
of both modalities as shift terms,

\begin{align}\label{eq:example}
\h(\ry \given \rx_1, \rx_2) = \hY(\ry) + \rx_1^\top\shiftparm_1 + \beta_2(\rx_2),
\end{align}

where $\hY:\calY \to \RR$ denotes the baseline transformation (the transformation
function obtained when $\rx_1^\top\shiftparm_1 + \beta_2(\rx_2) = 0$). Further,
$\shiftparm_1$ denotes the coefficients of the linear predictor and $\beta_2 :
\calX_2 \to \RR$ denotes the unstructured predictor for the text data, which are
typically controlled by a neural network. A  with (<ref>) is
distribution-free because for any constellation of covariates for which the shifting
predictor is zero, for all $(\rx^0_1, \rx^0_2) \in S^0 \coloneqq \{(\rx_1, \rx_2)
\in \calX_1 \times \calX_2 \given \rx_1^\top\shiftparm_1 + \beta_2(\rx_2) = 0\}$, and all
conditional distributions $Y \given \rX_1 = \rx^0_1, \rX_2 = \rx^0_2$, there
exists a unique baseline transformation given by $\hY = \pZ^{-1} \circ F_{\rY \given \rX_1 = \rx^0_1,
\rX_2 = \rx^0_2}$. In (<ref>), covariate effects are assumed to enter additively
on the scale of the transformation function, thus restricting distributions that
can be modeled for $(\rx_1, \rx_2) \in \calX \backslash S^0$. This argument can be extended to
more complex s <cit.>.
The example in (<ref>) is depicted in Figure <ref> for typical
types of responses and standard logistic latent distribution.

Example of a  with transformation function
$\h(\ry \given \rx_1, \rx_2)$ depending on a tabular modality $\rx_1$ and
a text modality $\rx_2$, which both enter as an additive shift term. The
tabular modality enters as a simple linear predictor $\rx_1^\top\shiftparm_1$
and the text data via the output of a neural network $\beta_2$, which is specified
by the user. Together with a baseline transformation $\hY$, whose parameterization
is discussed later, and the latent distribution $\pZ$, the DCTM is fully
specified. On the bottom, the transformation function $\h$ and cumulative distribution
function $F_{\rY \given \rX_1 = \rx_1, \rX_2 = \rx_2} = \pZ \circ \h$
are depicted for a binary, ordered, count, and continuous response for two realizations
of the tabular and text modalities. For binary and ordered responses with $K$ levels,
the transformation function contains one and $K - 1$ parameters, respectively, because
the CDF is constrained to one for the largest class.

In , the most general transformation function is parameterized
via basis expansions in the response and input modalities,

\begin{align} \label{eq:trafo}
\h(\ry \given \rx_{1:J}) = \left(\basisy(\ry) \otimes \basisx(\rx_{1:J})
\right)^\top\parm + \svec(\rx_{1:J})^\top\shiftparm, \quad \ry \in \calY,
\ \rx_{1:J} \in \calX,
\end{align}

where $\otimes$ denotes the Kronecker product and $\basisy, \basisx, \svec$ denote
the bases for the response, and the $J$ predictors, which either interact ($\basisx$)
with the response or simply shift ($\svec$) the transformation function. Thus,
shift effects are constant across all values of the response, the transformation
$\h$ can only shift up- or downwards (see Fig. <ref>). The effect
of interacting predictors may vary with the response and thus the shape of
$\h$ may change for different predictor values. For instance, an interacting
binary predictor leads to two separate transformations for each level, much like
stratum variables in survival analysis allow for separate hazard functions
[Collett, 2015].
However, in its general form, interacting predictors may also include neural
networks and thus unstructured predictors, making them extremely versatile.
Scale effects as introduced in Siegfried et al., 2022 are a special case of
interacting predictors, which are included in  by using $\basisx :
\calX \to \RR_+$ with $\rx_{1:J} \mapsto \sqrt{\exp(\gamma(\rx_{1:J}))}$ and, a neural network $\gamma : \calX \to \RR$. With a linear basis $\basisy$ in $\ry$,
$\basisy(\ry) = (1, \ry)^\top$, this is equivalent to location-scale regression with
error distribution $\pZ$.

Several types of univariate, potentially censored, responses can be handled. This
includes continuous ($\calY \subseteq \RR$), survival ($\calY \subseteq \RR_+$), count
($\calY = \mathbb{N}$), and ordered ($\calY = \{\ry_1, \dots, \ry_K\}$) responses.
For continuous responses, the basis of $\ry$ is a smooth function parameterized via
polynomials in Bernstein form of order $P$, denoted by $\bern{P}(\ry)$. For count
responses, the polynomials in Bernstein form are evaluated only at the integers,
$\bern{P}(\lfloor\ry\rfloor)$ [Siegfried and Hothorn, 2020]. For ordered responses, a
dummy-encoding is used, $\bvec(\ry_k) = \evec_k$, where $\evec_k$ denotes
the $k$-th unit vector. Linear and log-linear bases are supported as well. In
Appendix <ref>, we describe how the user can supply custom basis

The transformation function is required to be monotonically non-decreasing in $\ry$.
Hence, depending on the choice of basis, the parameters $\parm$ in (<ref>)
need to fulfill positivity or monotonicity constraints [Hothorn et al., 2014],
which can be enforced by appropriate reparameterizations. Without interacting predictors,
Bernstein polynomials and discrete bases require $\eparm_1 \leq \eparm_2 \leq \dots \leq
\eparm_M$ and linear and log-linear bases require positive slopes. For more complex
interacting predictors, the positivity of $\basisx$ has to be enforced together with
more complex constraints on $\parm$ to ensure a monotonically non-decreasing
transformation function <cit.>.

Finally, transformation models can be fitted by minimizing the negative average

\begin{align}\label{eq:nll}
\NLL \coloneqq - \frac{1}{n} \sum_{i = 1}^n \ell(\h; \ry_i, \rx_{1:J, i}),
\end{align}

where the observations $\{(\ry_i, \rx_{1:J,i})\}_{i=1}^n$ are assumed to be (conditionally)
independent. In , the default optimizer is (stochastic) gradient descent
using Adam [Kingma and Ba, 2015]. However, any keras [Allaire and Chollet, 2022] or
tensorflow optimizer
or a custom optimization routine can be used instead. For a single observation $(\ry,
\rx_{1:J})$, the log-likelihood contribution $\ell(\h; \ry, \rx_{1:J})$ depends on the
type of censoring of the observed response. Exact responses $\ry$ contribute
$\log\dZ(\h(\ry\given\rx))\h'(\ry \given \rx)$ to the log-likelihood. Interval-censored
responses $(\ubar\ry, \bar\ry]$ contribute $\log(\pZ(\h(\bar\ry\given\rx)) -
\pZ(\h(\ubar\ry \given \rx)))$. Left- and right-censored observations follow from the
interval-censored contribution as a special case, by letting $\ubar\ry \to -\infty$
and $\bar\ry \to +\infty$, respectively [Hothorn et al., 2014].

§.§ Autoregressive transformation models

Time series data pose one particular case where the independence assumption between observations
is not tenable and needs to be taken into
account. Formally, the joint distribution of a time series $(Y_t)_{t\in \mathcal{T}}$ with
$\mathcal{T} \subseteq \mathbb{N}_0$ can always be factorized in its conditional distributions,
by conditioning $Y_t$ on its full history $\calF_{t,1} \coloneqq (Y_{t-1}, \ldots, Y_{1})$.
A simplification is to impose a Markov property of order $p \geq 1$ which implies that the
conditional distribution of $Y_t$ only depends on the history up to and including $t-p$, that
is $\calF_{t,t-p} \coloneqq (Y_{t-1}, \ldots, Y_{t-p})$ rather than the entire history

Package  offers three ways on how to model time series data assuming the Markov
property. The naive way is given by classical transformation models where
$\calF_{t,t-p}$ is regarded in the basis expansion of the transformation function shown
in (<ref>) where elements of $\calF_{t,t-p}$ may interact with the response
$Y_t$ and simultaneously shift the transformation function. Furthermore,
[Rügamer et al., 2021] proposed the class of autoregressive transformation models (ATMs)
which differ from the naive approach (classical transformation models) in two perspectives.
First, the transformation function $h_t$ in ATMs can be time-varying which may result in different
transformations for different sub-periods. Second, the same $h_t$ is applied to $Y_t$ and each
element of $\calF_{t,t-p}$ simultaneously, resulting in a shared transformation between
$Y_t$ and its lags.

A special subclass of ATMs are AT($p$) models which do not allow for
interacting elements of $\calF_{t,t-p}$ with $Y_t$ through $\basisx$ but restrict to a
linear shift impact of the transformed
values of $\calF_{t,t-p}$ on the scale of $h$. The class of AT($p$) models is closely
related to a well-known class of time series models, autoregressive models of order $p$
<cit.>. In fact, AT($p$) models are equivalent to AR($p$) models
for $P=1$, $\svec(\rx_{1:J}) \equiv \rx_{1:J}$ and the independent white noise follows
the distribution $F_Z$ in a location-scale family <cit.>.
Learning the transformation simultaneously for the response and its lags as it is done in
ATMs is particularly important for ordinal time series, for which the dimensionality of the
model can thereby be reduced. Instead of modeling each level of the lagged response, the
one-dimensional transformed lagged response is included.
It also allows for a more consistent interpretation in the sense of autoregression because
we model $h(Y_t)$ at the current step (auto)regress the next time point  $h(Y_{t+1})$ on
the likewise transformed response $\h(\rY_t)$, not on the untransformed $Y_t$.
We showcase the practical differences between linear transformation models, AT($p$)
and ATM models in Section <ref>.

§.§ Application datasets

Movies data
We will illustrate the features of  on the movies dataset
[Kaggle, 2017]. The dataset contains information on 45,000 movies released prior
to July 2017, including number of ratings, budget, revenue, popularity, run time,
and genre. In addition, non-tabular reviews of the movies are available as text
data. We pre-process budget, revenue, and popularity using $\log(1 + x)$, due to
their skewed nature. For the text data, we use a text_tokenizer with a
1,000 word vocabulary, convert text to sequence and pad sequences to a maximum
length of 100 and truncate the end of a review.

Temperature data An application of autoregressive transformation
models to a time series of monthly mean maximum temperature in Melbourne (Australia)
in degrees Celsius between January 1971 and December 1990 (240 records) is
presented in Section <ref>. The temperature time series was recorded
by the Australian Bureau of Meteorology and later provided in Hyndman and Yang, 2022.

§ THE PACKAGE

Package  builds upon tensorflow as a fitting engine and
deepregression for setting up deep conditional transformation models.
In , response, interacting, and shifting terms are represented
as formula objects and correspond to the bases in (<ref>).
Internally, a processor is defined for each model term, which evaluates
its basis functions and optional penalties. For instance, for a continuous
response, a polynomial basis in Bernstein form and its derivatives are set up
by default (Table <ref>). Package  can include
terms modeled by user-specified neural network architectures for the
interacting and shifting terms (see Figure <ref> and
Figure <ref>). In the end, a single end-to-end trainable neural network
is set up, which may contain different neural network components for different
terms in the interacting or shifting predictor. Together with the supplied
latent_distr $\pZ$, the DCTM is fully specified and its parameters can be
estimated by minimizing the NLL via stochastic gradient descent (SGD). An
appropriate last-layer transformation ensures monotonicity constraints of the
interacting model term in the response. Each step in the  workflow
is highly customizable, custom functions for basis evaluation
(Appendix <ref>), custom last-layer transformations, and general-purpose
optimization routines (Section <ref>), such as SGD with adaptive learning
rates, can be supplied.

Schematic depiction of setting up and fitting DCTMs. Bases for
input predictors $\rx_{1:J}$ and response $\ry$ (circles) are evaluated and
enter the two neural network components interacting and
shifting according to (<ref>). The components'
outputs make up the transformation function $\h$. Together with the latent
distribution $\pZ$, the loss, NLL, can be evaluated and neural network
parameters can be optimized. Since $\pZ$ is parameter-free, all trainable
parameters are in the transformation function, as indicated by the dotted box.

Model function     Model name     Default basis     Default latent distribution

BoxCoxNN     Transformed normal     Bernstein     Standard normal

ColrNN     Continuous outcome logistic     Bernstein     Standard logistic

cotramNN     Count transformation     Bernstein     Standard logistic

CoxphNN     Cox proportional hazards     Bernstein     Standard minimum extreme value

LehmannNN     Lehmann-type     Bernstein     Standard maximum extreme value

LmNN     Normal linear     Linear     Standard normal

PolrNN     Proportional odds logistic     Discrete     Standard logistic

SurvregNN     Weibull     Log-linear     Standard minimum extreme value

Supported models together with the default choice of basis function
and latent distribution. Model functions summarized here are implemented
with a specific choice of basis and latent distribution that define
commonly applied regression models.

§.§ Main components

We describe the main components of  below by showing how to use
the formula interface, set up a DCTM, and fit the model. In this section, all
steps are illustrated with the movies dataset.

§.§.§ Formula interface

Models can be specified via a formula interface akin to the one used in
tram [Hothorn et al., 2022], where covariates interacting with the response
are supplied on the left-hand side, and shift effects are supplied on
the right-hand side of the formula, as illustrated below.

response | interacting   shifting

Thus, the formula interface mimics the transformation function as introduced in

\begin{align}\label{eq:annotatedtrafo}
\big(
\underbrace{\basisy(\ry)}_{\mathtt{response}}
\otimes
\overbrace{\basisx(\rx_{1:J})}^{\mathtt{interacting}}
\big)^\top\parm +
\end{align}

Case study: Formula interface

We begin by modeling the conditional distribution of vote_count
given a binary indicator of whether the movie is categorized as an action
movie or not (genreAction), the movie's budget and its
popularity. The below formula allows for separate baseline
transformations of the response for action movies vs. all other
genres, a smooth effect for budget and a linear effect for
popularity. Here, we use the standard spline basis representation
implemented in mgcv [Wood, 2021]. A non-exhaustive list of smoothers and other
processors is given in Table <ref>. Processors are
specialized functions for handling predictors which can speed up
computation. For instance, fac_processor() from safareg
evaluates factors on-line and thus facilitates computation for large
factor models <cit.>.

R> fm <- vote_count | genreAction   0 + s(budget, df = 3) + popularity

In the above formula we exclude an additional intercept in the shift
term by specifying 0 + ..., because the interacting basis
is already an intercept function.

Effect / Processor      Example formula
Linear      y x

Smooth      y s(x, ...)

Tensor product splines      y [te|ti|t2](x, ...)

Lasso      y lasso(x)

Group lasso      y grlasso(x)

Row-wise tensor product      y rwt(x)

Neural network      y nn(x)

Processor      *_processor

Implemented choices of interacting and shift
processors. All splines from mgcv are supported. Custom neural networks
can be supplied as functions via . Additional processors, for example,
for faster processing of large factors or interactions from safareg, can
be included via additional_processors
[Rügamer et al., 2022, Rügamer, 2022]. All terms can also be
included as interacting effects on the left-hand side of the formula, y | term(x, ...) 1.

§.§.§ Setting up DCTMs

DCTMs can be generically set up using the deeptrafo() function.

deeptrafo(formula = response | interacting   shifting, data = ...)

The data can be supplied as a data.frame or list. The
function returns a deeptrafo object, whose methods are described in
Section <ref>.

Special cases of DCTMs coincide with well-known models and are given their
own function in . The naming conventions in
follow the tram package [Hothorn et al., 2022] and add the suffix NN.
For instance, the proportional odds logistic regression model (ordinal response
and $\pZ = \expit$) is implemented as Polr() in tram and PolrNN()
in  (see Table <ref> for an overview).

Case study: Setting up DCTMs

For the movies data, we set up a count transformation model with standard
logistic latent distribution. We supply a custom Adam optimizer [Kingma and Ba, 2015]
for SGD with a learning rate of 0.1 decaying with a rate of $4 \cdot 10^{-4}$. The
training data train is the result of preprocessing steps described in
Section <ref>. The code for reproducing all output and figures can be
found on GitHub at <https://github.com/LucasKook/case-study-deeptrafo.git>.

R> opt <- optimizer_adam(learning_rate = 0.1, decay = 4e-4)
R> (m_fm <- cotramNN(formula = fm, data = train, optimizer = opt))

Count outcome deep conditional transformation model

Interacting:  vote_count | genreAction

Shifting:   0 + s(budget, df = 6) + popularity

Shift coefficients:
s(budget, df = 6)1 s(budget, df = 6)2 s(budget, df = 6)3 s(budget, df = 6)4
0.557             -0.702              0.760             -0.181
s(budget, df = 6)5 s(budget, df = 6)6 s(budget, df = 6)7 s(budget, df = 6)8
-0.201             -0.687              0.670              0.671
s(budget, df = 6)9         popularity
-0.377             -0.888

Printing the model to the console shows the model specification and shift
coefficients. Note that the model has only been initialized and not yet fitted.

§.§.§ Fitting DCTMs

For fitting DCTMs the user calls fit(), which calls the model
internal mod$fit_fun(), per default a wrapper around
fit.keras.engine.training.Model(), with the supplied arguments
(for instance epochs, batch_size). All functionalities
of fitting keras models carry over to fitting s, including
callbacks (custom operations applied after every iteration or mini-batch

Case study: Fitting DCTMs

The  object returned by cotramNN is fitted
for 1,000 epochs, with a batch size of 64, and a 10% validation split.
Below, we print the now-fitted model.

R> m_fm_hist <- fit(m_fm, epochs = 1e3, validation_split = 0.1, batch_size = 64,
+    verbose = FALSE)
R> unlist(coef(m_fm, which = "shifting"))

s(budget, df = 6)1 s(budget, df = 6)2 s(budget, df = 6)3 s(budget, df = 6)4
0.38339           -0.28824           -0.04608           -0.03992
s(budget, df = 6)5 s(budget, df = 6)6 s(budget, df = 6)7 s(budget, df = 6)8
0.00616           -0.02692           -0.00511            0.01355
s(budget, df = 6)9         popularity
-0.36587           -0.82771

Figure <ref>A depicts the training and validation loss trajectory
for inspecting convergence and overfitting, which can be generated with
plot(m_fm_hist). Figure <ref>B shows the estimated
transformation function. In Section <ref>, we describe how to produce
plots of the transformation function and density. Since genreAction is
included as a response-varying effect, the two transformation functions are
allowed to cross.

A: Training and validation loss trajectory for m_fm.
B: Estimated transformation functions for both levels of
genreAction with popularity and
budget fixed at their mean in the training

§.§.§ Working with neural networks

In s, neural networks map from a complex input space, such as text or
images, to a Euclidean space. When the neural network enters as a shift term,
the output of the network is a single real number which is interpretable on the
latent scale $\pZ^{-1}$, the scale of the transformation function. Custom
neural networks can be supplied to deeptrafo as functions or
keras_models via the  argument.

Case study: Working with neural networks

In our running example, we use the following architecture to model the contribution of
the movie reviews provided as textual descriptions. In Section <ref>, we
present an application with further downstream analysis of the text embedding.

R> embd_mod <- function(x) x |>
+    layer_embedding(input_dim = nr_words, output_dim = embedding_size) |>
+    layer_lstm(units = 50, return_sequences = TRUE) |>
+    layer_lstm(units = 50, return_sequences = FALSE) |>
+    layer_dropout(rate = 0.1) |> layer_dense(25) |>
+    layer_dropout(rate = 0.2) |> layer_dense(5) |>
+    layer_dropout(rate = 0.3) |> layer_dense(1)

The neural network embd_mod maps movie ratings to a real value (for more
details see the case study in Section <ref>). The interpretational
scale of output depends on the choice of latent distribution. Here, the
logistic distribution ($\pZ = \expit$) renders the output of embd_mod
interpretable on the log-odds scale. In turn, differences in the output of
embd_mod can be interpreted as log odds-ratios when changing, for
instance, a single word in a sentence and leaving everything else constant.
In our  model, we can now supply a named list
list(deep = embd_mod) and use deep(texts) in the formula.

R> fm_deep <- update(fm, .   . + deep(texts))
R> m_deep <- cotramNN(fm_deep, data = train,
+    list_of_deep_models = list(deep = embd_mod))
R> fit(m_deep, epochs = 50, validation_split = 0.1, batch_size = 32,
+    callbacks = list(callback_early_stopping(patience = 5)), verbose = FALSE)

Note that the default optimization routine may not produce optimization paths as
smooth as when omitting the neural network component. However, adaptively scheduled learning
rates for SGD often work well out-of-the-box, using optimizer = optimizer_adam()
as an argument when initializing the  model. Sometimes also different
learning schedules are needed for the different modalities (see Section <ref>).

§.§.§ Ensembling DCTMs

A simple and popular method to improve prediction performance and enable
uncertainty quantification are deep ensembles [Lakshminarayanan et al., 2017].
In a deep ensemble, a neural network model is trained $B$ times using the same training
and validation data, but different initial weights. Training via SGD may then converge
to different (local) minima and the members may yield different predictions. However,
averaging the predicted densities of the $B$ ensemble members is guaranteed to improve
upon the average individual performance (in terms of NLL). In ,
an ensemble of a model can be fitted via ensemble(). Besides classical deep
ensembling,  implements transformation ensembles [19].
Transformation ensembles are specifically tailored towards DCTMs and preserve their
additive structure and thus (partial) interpretability by averaging the predicted
transformation functions instead of the predicted densities.

Case study: Ensembling s

Epistemic uncertainty in the estimated smooth partial effect
of budget on vote_count.

Below, we fit five instances of m_deep. Then, we combine their
predictions on the scale of the transformation function and can investigate
uncertainty in the effects of the shifting predictors and prediction performance
on the test set.

R> ens_deep <- ensemble(m_deep, n_ensemble = 3, epochs = 50, batch_size = 64,
+    verbose = FALSE)

Figure <ref> shows the estimated smooth effect of budget
with pointwise epistemic uncertainty. Investigating the out-of-sample prediction
performance, we see that the transformation ensemble performs better than the
members do on average.

R> unlist(logLik(ens_deep, convert_fun = mean))

members1 members2 members3     mean ensemble
8.72     8.61     8.61     8.65     8.53

§.§.§ Cross-validating DCTMs

With cv(),  provides a cross-validation function for
objects. When supplying an integer to cv_folds, the
data is split into cv_folds number of folds. Alternatively, the user
can specify a list with two elements indicating data indices for training
and validation. The output of cv() can be used for tuning
smoothing hyperparameters, choosing between including a predictor as
interacting or shifting, or different neural network architectures.

Case study: Cross-validating s

The following call to cv() performs 5-fold cross validation while
fitting each instance of m_deep for 50 epochs. Train and validation
loss trajectories are shown in Figure <ref>. The vertical bars
indicate the epoch with the best average train/validation loss.

R> cv_deep <- cv(m_deep, epochs = 50, cv_folds = 5, batch_size = 64)
R> plot_cv(cv_deep)

Default plot generated by cv.deeptrafo(). The vertical
lines indicate the epoch with minimal average train/validation loss.

§.§ Methods overview

In the following, we briefly describe S3 methods for deeptrafo
and dtEnsemble objects.

§.§.§ Methods for deeptrafo objects

* coef return coefficients for the interacting
or shifting terms (controllable via
which_param = c("shifting", "interacting", "autoregressive")).

* predict returns in-sample predictions when
newdata is not supplied. The supported types are "trafo",
"pdf", "cdf", "interaction", "shift",
"terms". When newdata is supplied, predictions are
evaluated at the response, if it is contained in newdata.
The response can be omitted from newdata to predict the whole
conditional distribution. Then, predictions are evaluated on a grid
of length K, which is automatically generated based on the
response's support in the training data set. A custom grid of response
values can be supplied via q, which overwrites K.

* logLik evaluates in- or out-of-sample
log-likelihoods. This can be useful for model criticism and
evaluating predictive performance, respectively. The argument
convert_fun controls how the individual NLL contributions
are summarized. The default is function(x) = -sum(x) to
compute the log-likelihood. Other common choices include identity
to obtain the individual NLL contributions, or mean to get the
average NLL.

* plot by default plots smooth components in
the shifting formula part. Data for plotting can be obtained
by setting only_data = TRUE. Smooth terms in
interacting can be plotted by setting
which_param = "interacting". In the same manner as in
predict, densities evaluated in-sample (type = "pdf"),
CDFs (or probability integral transforms, with type = "cdf"),
and transformation functions (type = "trafo") can be obtained.
When omitting the response from newdata, the whole density,
cumulative distribution, or transformation function can be plotted.

* print prints a brief summary of the
including coefficients of additive linear and smooth terms in
shifting. Setting with_baseline = TRUE also prints
coefficients of linear and smooth terms in interacting.
The print_model argument toggles whether the keras
summary of the  should be printed in addition.

§.§.§ Methods for dtEnsemble objects

Methods coef and predict of deeptrafo objects take the same
arguments as their deeptrafo counterparts. The output is returned for
all ensemble members. Likewise, logLik returns the processed NLL
contributions for individual ensemble members, their average, and the
transformation ensemble.

§ APPLICATION: BINARY CLASSIFICATION

In this application, we use the movies dataset and fit four different
models with the goal to predict the binary response action (0: non-action
movie, 1: action movie, defined in the next code chunk),
which encodes whether a movie is an action movie or not. The model m_0
is unconditional; m_tab uses only one tabular predictor,
popularity, as linear shift predictor; m_text uses only
texts as an unstructured shift predictor; m_semi is a
semi-structured model which uses both modalities as shift predictors. First,
we encode the binary response as an ordered factor allowing us to use the
framework of ordinal neural network transformation models
[Kook et al., 2022].

R> train$action <- ordered(train$genreAction)
R> test$action <- ordered(test$genreAction, levels = levels(train$action))
\end{CodeInput}
\end{CodeChunk}

We then set up the formulas for the four models. The unconditional model is
specified without any predictors and \code{1} on the right-hand side. Later,
we will restrict this additional intercept to zero for identification (see
\code{warmstart_weights} in the definition of \code{m_0}). For all other
models, we remove the intercept directly by specifying \code{0 + ...} on
the right-hand side.
\begin{CodeChunk}
\begin{CodeInput}
R> fm_0 <- action ~ 1
R> fm_tab <- action ~ 0 + popularity
R> fm_text <- action ~ 0 + deep(texts)
R> fm_semi <- action ~ 0 + popularity + deep(texts)
\end{CodeInput}
\end{CodeChunk}
Note that \code{deep} is the same neural network architecture with text
embedding as in Section~\ref{subsec:maincomp}. Here, we use a custom
\cls{keras\_model} to which we can refer to by the name \code{"embd"}.
\begin{CodeChunk}
\begin{CodeInput}
R> make_keras_model <- function() {
+    return(keras_model_sequential(name = "embd") |>
+    layer_embedding(input_dim = nr_words, output_dim = embedding_size) |>
+    layer_lstm(units = 50, return_sequences = TRUE) |>
+    layer_lstm(units = 50, return_sequences = FALSE) |>
+    layer_dropout(rate = 0.1) |> layer_dense(25) |>
+    layer_dropout(rate = 0.2) |>  layer_dense(5, name = "penultimate") |>
+    layer_dropout(rate = 0.3) |> layer_dense(1))
+  }
\end{CodeInput}
\end{CodeChunk}

Next, we use \code{PolrNN()} to set up the different models with a standard
logistic latent distribution. Models including text data are trained
for ten epochs with early stopping and a patience of two, and the
weights from the epoch with the best validation loss are restored. The
unconditional and tabular-only models are trained full-batch and without
validation split until converging to the minimum.

\subsection{Unconditional model}

For the unconditional model, the intercept is fixed to zero via
\code{warmstart_weights} to ensure identification. The details
explaining the next code chunk can be found in Appendix~\ref{app:warm}.
\begin{CodeChunk}
\begin{CodeInput}
R> m_0 <- PolrNN(fm_0, data = train, optimizer = optimizer_adam(
+    learning_rate = 1e-2, decay = 1e-4), weight_options = weight_control(
+    general_weight_options = list(trainable = FALSE, use_bias = FALSE),
+    warmstart_weights = list(list(), list(), list("1" = 0))))
R> fit(m_0, epochs = 3e3, validation_split = 0, batch_size = length(
+    train$action), verbose = FALSE)
\end{CodeInput}
\end{CodeChunk}

The unconditional model \code{m_0} has one parameter which estimates the log-odds
of a movie belonging to a non-action genre without any predictors. The estimated intercept
parameter, given by \code{coef(m_0, which = "interacting")}, corresponds to the single
(fix) value of the transformation function $\h$ for a binary response
(see Figure~\ref{fig:extrafo}). The code chunk below shows that
the estimated intercept agrees with the observed log-odds of a movie belonging
to a non-action genre up to numerical inaccuracies.
\begin{CodeChunk}
\begin{CodeInput}
R> all.equal(unlist(unname(coef(m_0, which = "interacting"))),
+    qlogis(mean(train$action == 0)), tol = 1e-6)
\end{CodeInput}
\begin{CodeOutput}
[1] TRUE
\end{CodeOutput}
\end{CodeChunk}

We can obtain the unconditional log-odds also using \code{predict()} with
\code{type = "trafo"}. From the estimated log-odds we can determine the
probability for a movie to belong to a non-action genre which matches
the prevalence of non-action movies in the train set. The prevalence of non-action
movies can also be computed directly by using the \code{predict} function and setting
the argument \code{type = "pdf"} and supplying \code{action = 0} in \code{newdata}.

\subsection{Tabular-only model}

Next, we set up and fit \code{m_tab} including \code{popularity} as a linear shift
\begin{CodeChunk}
\begin{CodeInput}
R> m_tab <- PolrNN(fm_tab, data = train, optimizer = optimizer_adam(
+    learning_rate = 0.1, decay = 1e-4))
R> fit(m_tab, epochs = 1e3, batch_size = length(train$action),
+    validation_split = 0, verbose = FALSE)
\end{CodeInput}
\end{CodeChunk}

We obtain the estimated linear shift parameter $\hat\beta$ of \code{m_tab} by
\code{coef(m_tab, which_param = "shifting")}. Here, the odds for
a movie to belong to genre action change by the factor $\exp(-\hat\beta)$
\begin{CodeChunk}
\begin{CodeInput}
R> exp(-unlist(coef(m_tab, which = "shifting")))
\end{CodeInput}
\begin{CodeOutput}
\end{CodeOutput}
\end{CodeChunk}
when the predictor \code{popularity} increases by one unit. Without flipping the
sign, the coefficient $\hat\beta$ represents a log-odds ratio for a movie
belonging to a non-action genre compared to genre action upon a one-unit change
in \code{popularity}. Thus, the interpretation of $\hat\beta$ depends on the
parameterization of the model, in particular, the sign of the shifting predictor.
In \pkg{\deeptrafo{}}, the shifting predictor is consistently parameterized with
a plus sign for all models, which may differ from other implementations of the
same model type (\eg generalized linear models in \pkg{MASS}, or TMs in \pkg{tram}).

\subsection{Text-only model}

We now define and fit \code{m_text} including only the tokenized movie reviews.
\begin{CodeChunk}
\begin{CodeInput}
R> embd <- make_keras_model()
R> m_text <- PolrNN(fm_text, data = train, list_of_deep_models = list(
+    deep = embd), optimizer = optimizer_adam(learning_rate = 1e-4))
R> fit(m_text, epochs = 10, callbacks = list(callback_early_stopping(
+    patience = 2, restore_best_weights = TRUE)), verbose = FALSE)
\end{CodeInput}
\end{CodeChunk}
Analogously to smooth partial effects, the differences between two shift estimates
resulting from two different text inputs can still be interpreted as log odds-ratios.

We now have a closer look at what the \code{embd_mod} has learned. The network
takes as input the words (encoded as indices). Here, we use a vocabulary (all
words in the data set) of 10000 words and limit each review text to a size of
100 words. Review texts which are shorter are prepended with zeros, longer
movie descriptions are cut after 100 words. All punctuation is removed.
The \code{layer_embedding} learns to embed the word indices into an
\code{embedding_size}-dimensional representation. The resulting word embeddings
of the text are the input sequence to an LSTM layer with a 50-dimensional
memory state. The second LSTM layer outputs the 50-dimensional state after the last
word in the text, which is then further processed by a fully connected neural
network including dropout to prevent overfitting.

We can now use the trained \code{embd} to extract and analyze the derived
latent features of the embedding of single words or whole texts. We can obtain
the embedding of a single word as the output of \code{layer_embedding()}.
If we use a whole review as input, the latent features in the layer
\code{"penultimate"} correspond to a five-dimensional representation of
the text embedding processed by subsequent layers.

\begin{figure}[!t]
\centering
\includegraphics[width=\textwidth]{embedding-pca.pdf}
\caption{%
The first two principal components of the \code{"embedding"} layer for
single words (left) and the lower-dimensional representation learned in
\code{"penultimate"} (right) for whole movie reviews in the test data.
On the left, the PCA is computed on the word embedding of the 1,000 most
frequent words (we display only the 100 most frequent, black dots), and
on the right based on the low-dimensional representation of the embedding
of the 888 full movie reviews contained in the test data.
\end{figure}

Figure~\ref{fig:figure_movie_pca} shows the first two components of a
principle component analysis (PCA) applied to the word embedding (left) and
to the features learned in the \code{penultimate} layer for whole reviews
(right). The left plot reveals, that words hinting at an action movie, have a
similar embedding, and are separated from words that are rather representative
of a romance movie. The plot on the right of Figure~\ref{fig:figure_movie_pca}
confirms that the features derived from the embedding are tailored to discriminate
action movies from other genres since latent features of action movies cluster
together and are fairly well separated from romantic movies.

\subsection{Semi-structured model}

Finally, we set up the most complex model \code{m_semi} which takes both data
modalities as input. To achieve efficient training of the tabular part and
avoid overfitting of the embedding network \code{emdb_semi} we use two different
learning rates for the structured and unstructured part of the model. Specifically,
we optimize the intercept (with name \code{"ia_1__2"}) and tabular shift predictor
(with name \code{"popularity_3"}) with a higher learning rate, than the embedding
model (\code{"embd"}). Note, that in the embedding model some layers are named
explicitly, the names for the other components can be obtained from the
\cls{keras\_model} summary by initializing and calling
\code{print(m_semi, print\_model = TRUE)}.

\begin{CodeChunk}
\begin{CodeInput}
R> embd_semi <- make_keras_model()
R> optimizer <- function(model) {
+    optimizers_and_layers <- list(
+      tuple(optimizer_adam(learning_rate = 1e-2),
+      get_layer(model, "ia_1__2")),
+      tuple(optimizer_adam(learning_rate = 1e-2),
+      get_layer(model, "popularity_3")),
+      tuple(optimizer_adam(learning_rate = 1e-4),
+      get_layer(model, "embd")))
+    multioptimizer(optimizers_and_layers)
+  }
R> m_semi <- PolrNN(fm_semi, data = train, list_of_deep_models = list(
+    deep = embd_semi), optimizer = optimizer)
R> fit(m_semi, epochs = 10, callbacks = list(callback_early_stopping(
+    patience = 2, restore_best_weights = TRUE)), verbose = FALSE)
\end{CodeInput}
\end{CodeChunk}

\subsection{Model comparison}

Comparing the prediction performance of the models (measured in terms of
NLL) indicates that mainly the text modality contains information for separating
action movies from other genres. However, for a
more reliable assessment of this statement, the training schedule should be
optimized further. We compute 95\% bootstrap confidence intervals as a simple
uncertainty measure for the test NLL.
\begin{CodeChunk}
\begin{CodeInput}
R> bci <- function(mod) {
+    lli <- logLik(mod, newdata = test, convert_fun = identity)
+    bt <- boot(lli, statistic = \(x, d) mean(x[d]), R = 1e4)
+    btci <- boot.ci(bt, conf = 0.95, type = "perc")$percent[1, 4:5]
+    c("nll" = mean(lli), "lwr" = btci[1], "upr" = btci[2])
+ }

R> mods <- list("unconditional" = m_0, "tabular only" = m_tab,
+   "text only" = m_text, "semi-structured" = m_semi)
R> do.call("cbind", lapply(mods, bci))
\end{CodeInput}
\begin{CodeOutput}
unconditional tabular only text only semi-structured
nll         0.531        0.516     0.437           0.423
lwr         0.501        0.486     0.390           0.372
upr         0.562        0.549     0.486           0.478
\end{CodeOutput}
\end{CodeChunk}
Like \code{m_tab} the model \code{m_semi} estimates a linear shift parameter for
\code{popularity} which can also be interpreted as a (conditional) log
odds-ratio. The parameter goes in the same direction as in \code{m_tab} but has
a reduced absolute value and is now intepretable as a conditional log-odds
ratio because the text information that is now additionally accounted for.
% To obtain (epistemic) uncertainty estimates for coefficients in semi-structured
% models, the ensemble approach introduced in Section~\ref{sec:package} can be
% applied.
\begin{CodeChunk}
\begin{CodeInput}
R> c("tabular only" = unlist(unname(coef(m_tab))),
+    "semi-structured" = unlist(unname(coef(m_semi))))
\end{CodeInput}
\begin{CodeOutput}
tabular only semi-structured
-0.43           -0.32
\end{CodeOutput}
\end{CodeChunk}
The presented case study is meant to showcase some functionality of the package
\pkg{\deeptrafo} for binary responses. Note, that a \code{PolrNN} model with an
ordinal response that has $K$ levels and yields $K-1$ values for a discrete
transformation function (see Figure~\ref{fig:extrafo}) can be interpreted
analogously, \eg linear shift terms are still interpreted as log odds-ratios
\citep[for details and more examples see][]{kook2020ordinal}.

\section{Application: Autoregressive transformation models}\label{subsec:ATM}

We now return to ATMs, first discussed in Section~\ref{sec:atm}. One special form
of ATMs are AT($p$) models. AT($p$) models assume a linear impact of the transformed
values of $\calF_{t,t-p}$ on the scale of $h$. Because the transformation
is the same as for the response, AT($p$) models thus learn a joint transformation
of the response and its lags. For an illustration of transformation models
applied to time series data, the \code{temperature} dataset is used. We aim to estimate
the conditional distribution of the monthly mean maximum temperature in degrees
Celsius (°C) in Melbourne (Australia) between January 1971 and December 1990. A
descriptive analysis of the time series shows a strong seasonal pattern. This
motivates the application of a flexible approach that allows modeling the quickly
changing moments of the conditional distribution across time.

In the following, we compare three different forms of autoregressive transformation
models. The most flexible model (ATM) includes the lags as interacting predictors
and transformed lags in the shift term. The AT(3) model only includes the transformed
lags in the shift term. Lastly, the naive \code{ColrNN} model (Colr) includes the lags
as an additive linear term. In addition, every model contains a shift effect for \code{month}.
We compare the three models based on their estimated transformation functions and conditional
densities. We start by creating a factor variable \code{month} for the calendar month as well
as the lags $Y_{t-p}$, $p = 1, 2, 3$ denoted by \code{y_lag_<p>}. AT($p$) lags are included
using the internal \code{atplag()} processor.
\begin{CodeChunk}
\begin{CodeInput}
R> lags <- c(paste0("y_lag_", 1:p, collapse = "+"))
R> atplags <- c(paste0("atplag(y_lag_", 1:p, ")", collapse = "+"))
\end{CodeInput}
\end{CodeChunk}
The formula for the ATM model is given as follows. We include all three lags as
interacting predictors on the left-hand side of the formula and specify the \code{atplag}s
on the right-hand side.
\begin{CodeChunk}
\begin{CodeInput}
R> (fm_atm <- as.formula(paste0("y |", lags, "~ 0 + month +", atplags)))
\end{CodeInput}
\begin{CodeOutput}
y | y_lag_1 + y_lag_2 + y_lag_3 ~ 0 + month + atplag(y_lag_1) +
atplag(y_lag_2) + atplag(y_lag_3)
\end{CodeOutput}
\end{CodeChunk}
For the AT(3) model, we include the transformed lags in the shift
but not in the interacting term.
\begin{CodeChunk}
\begin{CodeInput}
R> (fm_atp <- as.formula(paste0("y ~ 0 + month +", atplags)))
\end{CodeInput}
\begin{CodeOutput}
y ~ 0 + month + atplag(y_lag_1) + atplag(y_lag_2) + atplag(y_lag_3)
\end{CodeOutput}
\end{CodeChunk}
The third model (Colr) we compare is a \code{ColrNN} model which includes
the raw lags as an additive shift term.
\begin{CodeChunk}
\begin{CodeInput}
R> (fm_colr <- as.formula(paste0("y ~ 0 + month + ", lags)))
\end{CodeInput}
\begin{CodeOutput}
y ~ 0 + month + y_lag_1 + y_lag_2 + y_lag_3
\end{CodeOutput}
\end{CodeChunk}
After preprocessing, the \code{temperature} dataset is saved in \code{d_ts}. We fix
the support of the response to \code{min_supp = 10} and \code{max_supp = 30} and
specify Bernstein polynomials of order \code{P = 6}. We use \code{ColrNN()}
to specify all models. ATM and AT(3) include \code{atplag}s and the third model, Colr,
does not.
\begin{CodeChunk}
\begin{CodeInput}
R> mod_fun <- function(fm) ColrNN(fm, data = d_ts,
+    trafo_options = trafo_control(order_bsp = P,
+    support = c(min_supp, max_supp)), tf_seed = 1,
+    optimizer = optimizer_adam(learning_rate = 0.01))
R> mods <- lapply(list(fm_atm, fm_atp, fm_colr), mod_fun)
\end{CodeInput}
\end{CodeChunk}
After defining the models, we proceed with training all three models.
In addition, we include callbacks to reduce the learning rate upon
encountering a plateau in the training loss, to ensure convergence of
the optimization procedure.
\begin{CodeChunk}
\begin{CodeInput}
R> fit_fun <- function(m) m |> fit(epochs = ep, callbacks = list(
+   callback_early_stopping(patience = 20, monitor = "val_loss"),
+   callback_reduce_lr_on_plateau(patience = 5)), batch_size = nrow(d_ts),
+   verbose = FALSE)
R> lapply(mods, fit_fun)
\end{CodeInput}
\end{CodeChunk}

\begin{figure}[!t]
\centering
\includegraphics[width=\textwidth]{plot-atps-atms-colrNN-1.pdf}
\caption{Estimated conditional densities (top row) of monthly temperature
records between June 1983 and May 1984, based on the ATM model (left), the
AT(3) model (center) and the Colr model (right). The observed
values across this time span are depicted in red. The plots in the bottom
row show the corresponding estimated conditional transformation functions.
\end{figure}

We compare the in-sample log-likelihood for the three models for the subset of
data between June 1977 and May 1978 in \code{t_idx}.
\begin{CodeChunk}
\begin{CodeInput}
R> t_idx <- seq(as.Date("1977-06-01"), as.Date("1978-05-01"), by = "month")
R> ndl <- d_ts[d_ts[["time"]] %in% t_idx]
R> structure(unlist(lapply(mods, logLik, newdata = ndl)),
+    names = c("ATM", paste0("AT(", p, ")"), "Colr"))
\end{CodeInput}
\begin{CodeOutput}
ATM AT(3)  Colr
-19.5 -22.5 -20.1
\end{CodeOutput}
\end{CodeChunk}
The comparison shows that the Colr and the ATM model fit similarly well
compared to the slightly less favorable fit of the AT(3) model. A visual
inspection of the estimated conditional densities depicted in
Figure~\ref{fig:figure_ATM_trafo} shows similar results for all three
estimation methods. In summary, the ATM class may be favored over naive
TMs (Colr) in the time series domain for its autoregressive structural
assumption, \ie lags entering in a transformed way, identical to the
transformation of $y_t$ \citep[see][]{rugamer2021timeseries}.

\section{Conclusion} \label{sec:conclusion}

With \pkg{\deeptrafo}, we introduce the first \proglang{R}~package for fitting a
broad class of distributional regression models with a neural network back-end.
Package \pkg{\deeptrafo} combines the advantages of transformation models, \ie flexible
distribution-free, yet interpretable models for conditional distributions, with
the advantages of neural network-based machine learning, which scales well for
large or non-tabular datasets. The intuitive formula interface allows users
familiar with packages such as \pkg{stats}, \pkg{MASS}, \pkg{tram}, \pkg{survival},
\pkg{mgcv}, and others to easily adapt their workflow to neural networks and more
complex datasets out-of-the-box. Users can supply custom basis functions,
optimization routines, and neural network architectures to tailor the \dctm{}
fully to their application. We illustrate \pkg{\deeptrafo} with tabular and text,
as well as time series data with count, discrete, and continuous outcomes, which
are all handled in a unified way. We demonstrate how custom neural network
architectures and optimizers can be used, and how to tune, evaluate, and interpret
\dctm{s}.

\section*{Acknowledgements}
We thank Sandra Siegfried for her valuable comments on the manuscript. The research
of DR has been partially supported by the German Federal Ministry of Education and
Research (BMBF) under Grant No. 01IS18036A.
The research of LK and BS was supported by Novartis Research Foundation (FreeNovation~2019)
and by the Swiss National Science Foundation (grant no. S-86013-01-01 and S-42344-04-01).
The authors of this work take full responsibilities for its content.

\begin{thebibliography}{34}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}

[Allaire and Chollet, 2022]
Allaire J, Chollet F (2022).
\newblock \emph{\pkg{keras}: R Interface to 'Keras'}.
\newblock \proglang{R} package version 2.9.0,
\urlprefix\url{https://CRAN.R-project.org/package=keras}.

[Allaire and Tang, 2022]
Allaire J, Tang Y (2022).
\newblock \emph{\pkg{tensorflow}: R Interface to 'TensorFlow'}.
\newblock \proglang{R} package version 2.9.0,
\urlprefix\url{https://CRAN.R-project.org/package=tensorflow}.

[Baumann \emph{et~al.}, 2021]
Baumann PFM, Hothorn T, Rügamer D (2021).
\newblock \enquote{Deep Conditional Transformation Models.}
\newblock In \emph{Machine Learning and Knowledge Discovery in Databases.
Research Track}, pp. 3--18. Springer-Verlag.
\newblock \doi{10.1007/978-3-030-86523-8_1}.

[Collett, 2015]
Collett D (2015).
\newblock \emph{Modelling Survival Data in Medical Research}.
\newblock CRC press.
\newblock \doi{10.1201/b18041}.

[Fahrmeir \emph{et~al.}, 2013]
Fahrmeir L, Kneib T, Lang S, Marx B (2013).
\newblock \emph{Regression - Models, Methods and Applications}.
\newblock Springer-Verlag, Berlin.

[Goodfellow \emph{et~al.}, 2016]
Goodfellow I, Bengio Y, Courville A (2016).
\newblock \emph{Deep Learning}.
\newblock MIT press.

[Hamilton, 2020]
Hamilton JD (2020).
\newblock \emph{Time Series Analysis}.
\newblock Princeton university press.
\newblock \doi{10.23943/princeton/9780691164502.003.0005}.

[Hothorn, 2020]
Hothorn T (2020).
\newblock \enquote{Transformation Boosting Machines.}
\newblock \emph{Statistics and Computing}, \textbf{30}(1), 141--152.
\newblock \doi{10.1007/s11222-019-09870-4}.

[Hothorn, 2021]
Hothorn T (2021).
\newblock \emph{\pkg{trtf}: Transformation Trees and Forests}.
\newblock \proglang{R} package version 0.3-8,
\urlprefix\url{https://CRAN.R-project.org/package=trtf}.

[Hothorn, 2022]
Hothorn T (2022).
\newblock \emph{\pkg{mlt}: Most Likely Transformations}.
\newblock \proglang{R} package version 1.4-3,
\urlprefix\url{https://CRAN.R-project.org/package=mlt}.

[Hothorn \emph{et~al.}, 2022]
Hothorn T, Barbanti L, Siegfried S (2022).
\newblock \emph{\pkg{tram}: Transformation Models}.
\newblock \proglang{R} package version 0.7-2,
\urlprefix\url{https://CRAN.R-project.org/package=tram}.

[Hothorn \emph{et~al.}, 2014]
Hothorn T, Kneib T, Bühlmann P (2014).
\newblock \enquote{{Conditional Transformation Models}.}
\newblock \emph{Journal of the Royal Statistical Society~B: Statistical
Methodology}, \textbf{76}(1), 3--27.
\newblock \doi{10.1111/rssb.12017}.

[Hothorn \emph{et~al.}, 2018]
Hothorn T, Möst L, Bühlmann P (2018).
\newblock \enquote{{Most Likely Transformations}.}
\newblock \emph{Scandinavian Journal of Statistics}, \textbf{45}(1), 110--134.
\newblock \doi{10.1111/sjos.12291}.

[Hyndman and Yang, 2022]
Hyndman R, Yang Y (2022).
\newblock \emph{\pkg{tsdl}: Time Series Data Library}.
\newblock \proglang{R}~package version 0.1.0,
\urlprefix\url{https://finyang.github.io/tsdl/}.

[Kaggle}, 2017]
{Kaggle} (2017).
\newblock \enquote{The \code{Movies} Dataset.}
\newblock
\urlprefix\url{https://www.kaggle.com/datasets/rounakbanik/the-movies-dataset}.

[Kingma and Ba, 2015]
Kingma DP, Ba JL (2015).
\newblock \enquote{{Adam: A Method for Stochastic Optimization}.}
\newblock In \emph{3rd International Conference on Learning Representations,
ICLR 2015 - Conference Track Proceedings}. International Conference on
Learning Representations, ICLR.
\newblock \doi{10.48550/arxiv.1412.6980}.

[17]
Kook L ({2022}).
\newblock \enquote{{\pkg{tramvs}: Optimal Subset Selection in Transformation
\newblock \proglang{R} package version 0.0-3,
\urlprefix\url{https://CRAN.R-project.org/package=tramvs}.

[Kook \emph{et~al.}, 2022]
Kook L, Baumann PFM, R\"ugamer D (2022).
\newblock \emph{{\pkg{deeptrafo}: Fitting Deep Conditional Transformation
\newblock \proglang{R} package version 0.1-1,
\urlprefix\url{https://CRAN.R-project.org/package=deeptrafo}.

[19]
Kook L, Götschi A, Baumann PFM, Hothorn T, Sick B ({2022}).
\newblock \enquote{Deep Interpretable Ensembles.}
\newblock \emph{arXiv Preprint arXiv:2205.12729}.
\newblock \doi{10.48550/arxiv.2205.12729}.

[Kook \emph{et~al.}, 2022]
Kook L, Herzog L, Hothorn T, Dürr O, Sick B (2022).
\newblock \enquote{Deep and Interpretable Regression Models for Ordinal
\newblock \emph{Pattern Recognition}, \textbf{122}, 108263.
\newblock \doi{10.1016/j.patcog.2021.108263}.

[Kook and Hothorn, 2021]
Kook L, Hothorn T (2021).
\newblock \enquote{{Regularized Transformation Models: The \pkg{tramnet}
\newblock \emph{{The \proglang{R} Journal}}, \textbf{13}(1), 581--594.
\newblock \doi{10.32614/rj-2021-054}.

[Lakshminarayanan \emph{et~al.}, 2017]
Lakshminarayanan B, Pritzel A, Blundell C (2017).
\newblock \enquote{Simple and Scalable Predictive Uncertainty Estimation Using
Deep Ensembles.}
\newblock In I~Guyon, UV~Luxburg, S~Bengio, H~Wallach, R~Fergus,
S~Vishwanathan, R~Garnett (eds.), \emph{Advances in Neural Information
Processing Systems}, volume~30. Curran Associates, Inc.
\newblock
\urlprefix\url{https://proceedings.neurips.cc/paper/2017/file/9ef2ed4b7fd2c810847ffa5fa85bce38-Paper.pdf}.

[Rügamer, 2022]
Rügamer D (2022).
\newblock \enquote{Additive Higher-Order Factorization Machines.}
\newblock \emph{arXiv Preprint arXiv:2205.14515}.
\newblock \doi{10.48550/arxiv.2205.14515}.

[Rügamer \emph{et~al.}, 2021]
Rügamer D, Baumann PFM, Kneib T, Hothorn T (2021{\natexlab{a}}).
\newblock \enquote{Probabilistic Time Series Forecasts with Autoregressive
Transformation Models.}
\newblock \emph{arXiv Preprint}.
\newblock \doi{10.48550/arxiv.2110.08248}.
\newblock To appear in \emph{Statistics \& Computing}.

[Rügamer \emph{et~al.}, 2022]
Rügamer D, Bender A, Wiegrebe S, Racek D, Bischl B, Müller C, Stachl C
\newblock \enquote{Factorized Structured Regression for Large-Scale Varying
Coefficient Models.}
\newblock \emph{arXiv Preprint arXiv:2205.13080}.
\newblock \doi{10.48550/arxiv.2205.13080}.

[Rügamer \emph{et~al.}, 2021]
Rügamer D, Kolb C, Fritz C, Pfisterer F, Bischl B, Shen R, Bukas C, Thalmeier
D, Baumann P, Kook L, Klein N, \emph{et~al.} (2021{\natexlab{b}}).
\newblock \enquote{\pkg{deepregression}: a Flexible Neural Network Framework
for Semi-Structured Deep Distributional Regression.}
\newblock \emph{arXiv Preprint}.
\newblock \doi{10.48550/arxiv.2104.02705}.
\newblock To appear in \emph{Journal of Statistical Software}.

[Rügamer \emph{et~al.}, 2020]
Rügamer D, Kolb C, Klein N (2020).
\newblock \enquote{Semi-Structured Deep Distributional Regression: Combining
Structured Additive Models and Deep Learning.}
\newblock \emph{arXiv Preprint}.
\newblock \doi{10.48550/arxiv.2002.05777}.

[Sick \emph{et~al.}, 2021]
Sick B, Hothorn T, Dürr O (2021).
\newblock \enquote{Deep Transformation Models: Tackling Complex Regression
Problems with Neural Network Based Transformation Models.}
\newblock In \emph{2020 25th International Conference on Pattern Recognition
({ICPR})}. {IEEE}.
\newblock \doi{10.1109/icpr48806.2021.9413177}.

[Siegfried and Hothorn, 2020]
Siegfried S, Hothorn T (2020).
\newblock \enquote{Count Transformation Models.}
\newblock \emph{Methods in Ecology and Evolution}, \textbf{11}(7), 818--827.
\newblock \doi{10.1111/2041-210x.13383}.

[Siegfried \emph{et~al.}, 2022]
Siegfried S, Kook L, Hothorn T (2022).
\newblock \enquote{Distribution-Free Location-Scale Regression.}
\newblock \emph{arXiv Preprint arXiv:2208.05302}.
\newblock \doi{10.48550/arxiv.2208.05302}.

[Tamási and Hothorn, 2021]
Tamási B, Hothorn T (2021).
\newblock \enquote{\pkg{tramME}: Mixed-Effects Transformation Models Using
Template Model Builder.}
\newblock \emph{The R Journal}, \textbf{13}(2), 398--418.
\newblock \doi{10.32614/rj-2021-075}.

[Varadhan, 2022]
Varadhan R (2022).
\newblock \emph{\pkg{alabama}: Constrained Nonlinear Optimization}.
\newblock \proglang{R} package version 2022.4-1,
\urlprefix\url{https://CRAN.R-project.org/package=alabama}.

[Varadhan and Gilbert, 2019]
Varadhan R, Gilbert P (2019).
\newblock \emph{\pkg{BB}: Solving and Optimizing Large-Scale Nonlinear
\newblock \proglang{R} package version 2019.10-1,
\urlprefix\url{https://CRAN.R-project.org/package=BB}.

[Wood, 2021]
Wood S (2021).
\newblock \emph{\pkg{mgcv}: Mixed GAM Computation Vehicle with Automatic
Smoothness Estimation}.
\newblock \proglang{R} package version 1.8-38,
\urlprefix\url{https://CRAN.R-project.org/package=mgcv}.

\end{thebibliography}

\appendix

\section{Appendix}

In this appendix, we describe how \pkg{\deeptrafo} handles censored responses
(Appendix~\ref{app:cens}), how the user can warmstart and fix weights of
interacting and shifting terms (Appendix~\ref{app:warm}), and how to include
custom basis functions (Appendix~\ref{app:basis}). In addition, we describe an
alternative formula interface (Appendix~\ref{app:form}) and show how to use
\pkg{\deeptrafo} for large tabular datasets (Appendix~\ref{app:fac}).

\section{Handling censored responses}\label{app:cens}

Package \pkg{\deeptrafo} detects the type of response automatically. However, the user
may specify the type explicitly via \code{response\_type} in \code{deeptrafo()}
and all alias/wrapper functions. Allowed types of responses are continuous,
count, survival, ordered (including binary). Censored responses can be supplied
as \cls{Surv} objects. Internally, ordered and count responses are treated
as censored. For instance, the two observations \code{c(0L, 1L)} with
\code{response\_type = "count"} are internally represented as left- and
interval-censored, respectively.

\begin{CodeChunk}
\begin{CodeInput}
R> deeptrafo:::response(y = c(0L, 1L))
\end{CodeInput}
\begin{CodeOutput}
cleft exact cright cinterval
[1,]     1     0      0         0
[2,]     0     0      0         1
[1] "count"
\end{CodeOutput}
\end{CodeChunk}

% \todo{@LK: allow response object from mlt}

\section{Warmstarting and fixing weights}\label{app:warm}

Warmstarting and fixing weights may be important in numerical experiments, for
finetuning parts of the models, or transfer learning [Goodfellow \emph{et~al.}, 2016].
In \pkg{\deeptrafo}, the user can supply a \cls{keras\_model}, as returned,
for instance, by \linebreak\code{keras_model_sequential()}. When defining the
model, \pkg{keras} specific arguments for controlling weight initialization
can be used, as shown below.
\begin{CodeChunk}
\begin{CodeInput}
R> nn <- keras_model_sequential() |>
+   layer_dense(input_shape = 1L, units = 3L, activation = "relu",
+   use_bias = FALSE, kernel_initializer = initializer_constant(
+   value = 1))
R> unlist(get_weights(nn))
\end{CodeInput}
\begin{CodeOutput}
[1] 1 1 1
\end{CodeOutput}
\end{CodeChunk}

To warmstart or fix coefficients of the interacting or shifting part of
a \dctm{}, the \code{weight\_options} argument in \code{deeptrafo()} can
supplied with the output of \code{weight\_control()}, which, in addition
to others, takes the same arguments as the \pkg{keras} layers above.
\begin{CodeChunk}
\begin{CodeInput}
R> args(weight_control)
\end{CodeInput}
\begin{CodeOutput}
function (specific_weight_options = NULL, general_weight_options = list(
activation = NULL, use_bias = FALSE, trainable = TRUE,
kernel_initializer = "glorot_uniform", bias_initializer = "zeros",
kernel_regularizer = NULL, bias_regularizer = NULL,
activity_regularizer = NULL, kernel_constraint = NULL,
bias_constraint = NULL), warmstart_weights = NULL,
shared_layers = NULL)
\end{CodeOutput}
\end{CodeChunk}

Below, we warmstart the shift coefficient for a \code{PolrNN} model.
Here, \code{warmstart\_weights} takes a list with three components,
of which the first two control the weights of the interacting predictor
and the last the weights of the shift predictor. The weights can be
referred to by the name of the covariate, \ie \code{"temp" = 0}.
\begin{CodeChunk}
\begin{CodeInput}
R> data("wine", package = "ordinal")
R> mw <- deeptrafo(
+   response ~ 0 + temp,
+   data = wine, weight_options = weight_control(warmstart_weights = list(
+     list(), list(), list("temp" = 0))))
R> unlist(coef(mw))
\end{CodeInput}
\begin{CodeOutput}
\end{CodeOutput}
\end{CodeChunk}
The three lists correspond to the three formula components
\code{response}, \code{interacting}, and \code{shifting}.
Note that the list corresponding to the response is always
empty, since it does not contain any parameters. In case there
is no interacting predictor, the second list corresponds to
the parameters of the basis function of the response, \ie
the intercept function. In case there is no shift term, an
intercept is set up which can be referred to as \code{"1"}
and frozen as illustrated in the main text (Section~\ref{sec:ontram}).
In the example above, we warmstart weights of a component
in the shift term and supply two empty lists for the other

\section{Including custom basis functions}\label{app:basis}

Linear, log-linear, and Bernstein bases, as used by \pkg{\deeptrafo},
require (linear) inequality constraints on their parameters. Internally,
these constraints are handled in \code{trafo\_control()}, by supplying
an \pkg{keras} layer, which transforms the weights for the interacting
predictor appropriately. In \pkg{\deeptrafo}, the implemented bases
are \code{"bernstein"}, \code{"ordered"}, and \code{"shiftscale"}. The
former two require $\eparm_{jP + 1)} \leq \eparm_{jP + 2} \leq \dots
\leq \eparm_{jP + P}$, $l = 0, \dots, L - 1$ for $\basisx(\rx_{1:J}) \in
\RR^L$ and degree $P - 1$ Bernstein basis or ordered response with $P + 1$
levels. The shift-scale basis requires only $\eparm_1 > 0$ in $\ry \mapsto
\eparm_0 + \eparm_1\ry$.

The user can now supply custom basis functions as shown below. First,
the basis (\code{linear_basis}) and its derivative
(\code{linear_basis_prime}) are defined. Afterwards, the constraints
on the parameters are defined using \proglang{Python}- and
\pkg{tensorflow}-specified constructs (\code{tf$...}).

\begin{CodeChunk}
\begin{CodeInput}
R> linear_basis <- function(y) {
+    ret <- cbind(1, y)
+    if (NROW(ret) == 1)
+      return(as.vector(ret))
+    ret
+  }
R> linear_basis_prime <- function(y) {
+    ret <- cbind(0, rep(1, length(y)))
+    if (NROW(ret) == 1)
+      return(as.vector(ret))
+    ret
+  }
R> constraint <- function(w, bsp_dim) {
+    w_res <- tf$reshape(w, shape = list(bsp_dim, as.integer(nrow(w) /
+    bsp_dim)))
+    w1 <- tf$slice(w_res, c(0L, 0L), size = c(1L, ncol(w_res)))
+    wrest <- tf$math$softplus(tf$slice(w_res, c(1L, 0L), size = c(
+    as.integer(nrow(w_res) - 1), ncol(w_res))))
+    w_w_cons <- k_concatenate(list(w1, wrest), axis = 1L)
+    tf$reshape(w_w_cons, shape = list(nrow(w), 1L))
+  }
R> tfc <- trafo_control(
+    order_bsp = 1L,
+    y_basis_fun = linear_basis,
+    y_basis_fun_prime = linear_basis_prime,
+    basis = constraint
+  )
\end{CodeInput}
\end{CodeChunk}

We can now compare our re-implementation of a transformation model with
linear basis against \code{Lm()} from \pkg{tram}. To efficiently fit
\dctm{s} for small tabular datasets, we recommend full-batch (\ie batch size
$n$) training with a large learning rate (0.01) and either decay or callbacks
for reducing the learning on validation loss plateaus.
\begin{CodeChunk}
\begin{CodeInput}
R> library("tram")
R> set.seed(1)
R> n <- 1e3
R> d <- data.frame(y = 1 + rnorm(1e3), x = rnorm(1e3))
R> m <- deeptrafo(y ~ 0 + x, data = d, trafo_options = tfc,
+    optimizer = optimizer_adam(learning_rate = 1e-2),
+    latent_distr = "normal")
R> fit(m, batch_size = 1e3, epochs = 3e3, validation_split = NULL,
+    callbacks = list(callback_reduce_lr_on_plateau(monitor = "loss")),
+    verbose = FALSE)
R> abs(unlist(coef(m)) - coef(Lm(y ~ x, data = d)))
\end{CodeInput}
\begin{CodeOutput}
\end{CodeOutput}
\end{CodeChunk}

\section{Alternative formula interface}\label{app:form}

Following {\sc ontram}s (ordinal neural network transformation models),
introduced in Kook \& Herzog \emph{et al.} ([Kook \emph{et~al.}, 2022]),
\pkg{\deeptrafo} offers an alternative formula interface. Here, the user
supplies a separate formula for the intercepts (before: interacting) and
for the shift (before: shifting) and avoids using the pipe \code{|} on
the left-hand-side of
the formula. Internally, the formula is translated back into the form in
(\ref{eq:trafo}). All other functionalities in the article carry over
to {\sc ontram}s as well. The same interface for other than ordinal
responses is implemented in \code{dctm()}.
\begin{CodeChunk}
\begin{CodeInput}
R> dord <- data.frame(Y = ordered(sample.int(6, 100, TRUE)),
+    X = rnorm(100), Z = rnorm(100))
R> ontram(response = ~ Y, intercept = ~ X, shift = ~ 0 + s(Z, df = 3),
+    data = dord)
\end{CodeInput}
\begin{CodeOutput}
Ordinal outcome deep conditional transformation model

Interacting:  Y | X

Shifting:  ~0 + s(Z, df = 3)

Shift coefficients:
s(Z, df = 3)1 s(Z, df = 3)2 s(Z, df = 3)3 s(Z, df = 3)4 s(Z, df = 3)5
-0.4760       -0.7326       -0.6233       -0.4061       -0.4309
s(Z, df = 3)6 s(Z, df = 3)7 s(Z, df = 3)8 s(Z, df = 3)9
-0.5447        0.6729        0.7376        0.0947
\end{CodeOutput}
\end{CodeChunk}

\section{Large factor models}\label{app:fac}

We consider a large factor model with $10^6$ observations and a factor
variable with $10^3$ levels. The standard implementation of \code{lm()}
and \code{LmNN()} fail to process the data, due to evaluating the large
model matrix. However, we can use \code{fac_processor()} from \pkg{safareg}
to circumvent this issue and use mini-batch stochastic gradient descent
to fit the model on a standard machine. Now, \pkg{\deeptrafo} can fit
large factor models for arbitrary types of responses and censoring.
\begin{CodeChunk}
\begin{CodeInput}
R> set.seed(0)
R> library("safareg")
R> n <- 1e6
R> nlevs <- 1e3
R> X <- factor(sample.int(nlevs, n, TRUE))
R> Y <- (X == 2) - (X == 3) + rnorm(n)
R> d <- data.frame(Y = Y, X = X)
R> m <- LmNN(Y ~ 0 + fac(X), data = d, additional_processor = list(
+    fac = fac_processor))
R> fit(m, batch_size = 1e4, epochs = 20, validation_split = 0,
+    callbacks = list(callback_early_stopping("loss", patience = 3),
+    callback_reduce_lr_on_plateau("loss", 0.9, 2)))
R> bl <- unlist(coef(m, which = "interacting"))
R> - (unlist(coef(m))[1:5] + bl[1]) / bl[2]
\end{CodeInput}
\begin{CodeOutput}
fac(X)1 fac(X)2 fac(X)3 fac(X)4 fac(X)5
-0.0204  0.9986 -1.0156 -0.0249  0.0477
\end{CodeOutput}
\end{CodeChunk}

\end{document}