Regularized Transformation Models: The tramnet Package

The tramnet package implements regularized linear transformation models by combining the ﬂexible class of transformation models from tram with constrained convex optimization implemented in CVXR . Regularized transformation models unify many existing and novel regularized regression models under one theoretical and computational framework. Regularization strategies implemented for transformation models in tramnet include the LASSO, ridge regression and the elastic net and follow the parametrization in glmnet . Several functionalities for optimizing the hyperparameters, including model-based optimization based on the mlrMBO package, are implemented. A multitude of S3 methods are deployed for visualization, handling and simulation purposes. This work aims at illustrating all facets of tramnet in realistic settings and comparing regularized transformation models with existing implementations of similar models.


Introduction
A plethora of R packages exist to estimate generalized linear regression models via penalized maximum likelihood, such as penalized (Goeman, 2010) and glmnet (Friedman et al., 2010). Both packages come with an extension to fit a penalized form of the Cox proportional hazard model. The tramnet package aims at unifying the above mentioned and several novel models using the theoretical and computational framework of transformation models. Novel models in this class include Continuous Outcome Logistic Regression (COLR) as introduced by Lohse et al. (2017) and Box-Cox type regression models with a transformed conditionally normal response (Box and Cox, 1964;Hothorn, 2020c).
The disciplined convex optimization package CVXR (Fu et al., 2020) is applied to solve the constrained convex optimization problems that arise when fitting regularized transformation models. Transformation models are introduced in Section 2.1.1, for a more theoretical treatise we refer to Hothorn et al. (2014Hothorn et al. ( , 2018; Hothorn (2020b). Convex optimization and domain specific languages are briefly discussed in Section 2.1.3, followed by a treatment of model-based optimization for hyperparameter tuning (2.1.4).

Transformation models
In stark contrast to penalized generalized linear models, regularized transformation models aim at estimating the response's whole conditional distribution instead of focusing on a single moment, e.g. the conditional mean. This conditional distribution function of a response Y is decomposed into an a priori chosen absolute continuous and log-concave error distribution F and a conditional transformation function h(y|x, s) that depends on the measured covariates x and stratum variables s and is monotone increasing in y. Although the model class is more flexible, packages tram and tramnet focus on stratified linear transformation models of the form P (Y ≤ y|X = x, S = s) = F (h(y|s, x)) = F h(y|s) − x β . (1) Here, the baseline transformation is allowed to vary with stratum variables s, while covariate effects β are restricted to be shifts common to all baseline transformations h(y|s).
In order for the model to represent a valid cumulative distribution function, F (h(y|s, x)) has to be monotone increasing in y and thus in h for all possible strata s and all possible configurations of the covariates x. To ensure monotonicity, h is parametrized in terms of a basis expansion using Bernstein polynomials as implemented in the basefun package (Hothorn, 2020b). Hence, h is of the form h(y) = a Bs,p (y) ϑ, where a Bs,p (y) denotes the vector of basis functions in y of order p and ϑ are the coefficients for each basis function. Conveniently, a Bs,p (y) ϑ is monotone increasing in y as long as holds. For the concrete parameterization of stratified linear transformation models the reader is referred to Hothorn (2020c).
Many contemporary models can be understood as linear transformation models, such as the normal linear regression model, logistic regression for binary, ordered and continuous responses, as well as exponential, Weibull and Rayleigh regression and the Cox model in survival analysis. Thus, by appropriately choosing and parametrizing F and h one can understand all those models in the same maximum likelihood-based framework. One can formulate the corresponding likelihood contributions not only for exact observations, but under any form of random censoring and truncation for continuous and discrete or ordered categorical responses.
Given a univariate response Y and a set of covariates X one can specify the following cumulative distribution function and density valid for any linear transformation model, From here, the log-likelihood contributions for exact, right, left, and interval censored responses can be derived as The joint log-likelihood of several observations y 1 , . . . , y n is obtained by summing over the individual log-likelihood contributions i under the assumption that the individual samples are independent and identically distributed, the case exclusively dealt with by tramnet.

Regularization
The aim of tramnet is to enable the estimation of regularized stratified linear transformation models. This is achieved by optimizing a penalized form of the log-likelihood introduced in the last section. The penalized log-likelihood, consists of the unpenalized log-likelihood and an additional penalty term. Note that only the shift parameters β are penalized, whereas the coefficients for the baseline transformation ϑ remain unpenalized. The parameterization of the penalty is chosen to be the same as in glmnet, consisting of a global penalization parameter λ and a mixing parameter α controlling the amount of L 1 compared to L 2 penalization.
The two penalties and any combination thereof have unique properties and may be useful under different circumstances. A pure L 1 penalty was first introduced by Tibshirani (1996) in an OLS framework and was dubbed the LASSO (Least Absolute Shrinkage and Selection Operator) due to its property of shrinking regression coefficients exactly to 0 for large enough λ. A pure LASSO penalty can be obtained in a regularized transformation model by specifying α = 1. Applying an L 2 penalty in an OLS problem was introduced more than five decades earlier by Tikhonov (1943) and later termed ridge regression (Hoerl and Kennard, 1970). In contrast to LASSO, ridge regression leads to shrunken regression coefficients, but does not perform automatic variable selection. Zou and Hastie (2005) picked up on both approaches, discussed their advantages, disadvantages and overall characteristics and combined them into the elastic net penalty, a convex combination of an L 1 and L 2 penatly controlled by the mixing parameter α. Some of these properties will be illustrated for different models and a real world data set in sections 2.1.6 and 2.2.2.

Constrained convex optimization
Special algorithms were developed to optimize regularized objective functions, most prominently the LARS and LARS-EN algorithm (Efron et al., 2004) and variants thereof for the penalized Cox model (Goeman, 2010). However, the aim of tramnet is to solve the objective functions arising in regularized transformation models in a single computational framework. Due to the log-concavity of all choices for F in this package and h(y) being monotone increasing in y, the resulting log-likelihood contributions for any form of censoring and truncation are concave and can thus be solved by constrained convex optimization.
The fairly recent development of CVXR allows the specification of constrained convex optimization problems in terms of domain-specific language, yielding an intuitive and highly flexible framework for constrained optimization. Because checking convexity of an aribitrarily complex expression is extremely hard, CVXR makes use of a library of smaller expressions, called atoms, with known monotonicity and curvature and tries to decompose the objective at hand using a set of rules from disciplined convex optimization (DCP) (Grant et al., 2006). Thus a complex expression's curvature can be more easily determined.
More formally, convex optimization aims at solving a problem of the form where ϑ ∈ R p is the parameter vector, g(ϑ) is the objective function to be optimized, g i (ϑ) specify the inequality constraints and A ∈ R n×p and b ∈ R p parametrize any equality constraints on ϑ. Importantly, the objective function and all inequality constraint functions are convex (Boyd and Vandenberghe, 2004).
The likelihood ∑ i i (ϑ, β; y i , s i , x i ) for transformation models of the form (1) are convex for error distributions with log-concave density, because log-convexity of −F ensures the existence and uniqueness of the most likely transformationĥ(y) and the convexity of − (h; y, x). Because the penalty term is convex in β, it can be added to the negative log-likelihood while conserving convexity. However, monotonicity of h imposes inequality constraints on the parameters of the baseline transformation as illustrated in equation (2). The elegance of domain-specific language based optimizers comes to play when adding these and potential other inequality or equality constraints to the objective function, which will be showcased in Section 2.2.3. Thus, the optimisation routines implemented in package CVXR can be applied for computing maximum likelihood estimates of the parameters of model (1).

Model-based optimization
The predictive capabilities of regularized regression models heavily depend on the hyperparameters α and λ. Hyperparameter tuning can be addressed by a multitude of methods with varying computational complexity, advantages and disadvantages. Naive or random grid search for more than one tuning parameter are computationally demanding, especially if the objective function is expensive to evaluate. Model-based optimization circumvents this issue by fitting a surrogate model, usually a Gaussian process, to the objective function. The objective function is evaluated at an initial, e.g. a random latin hypercube, design, to which the Gaussian process is subsequently fit. The surrogate model then proposes the next set of hyperparameters to evaluate the objective function at by some infill criterion (Horn and Bischl, 2016). Bischl et al. (2017) implement model-based optimization for multi-objective blackbox functions in the mlrMBO package. The objective function can in theory be vector-valued and the tuning parameter spaces may be categorical. In tramnet the objective function is the cross-validated log-likelihood optimized using a Kriging surrogate model with expected improvement as the infill criterion. Model-based optimization for hyperparameter tuning is illustrated in section 2.2.

Basic usage
The initial step is fitting a potentially stratified, transformation model of the form R> m1 <-tram(y | s~1, ...) omitting all explanatory variables. This sets up the basis expansion for the transformation function, whose regression coefficients will not be penalized, as mentioned in section 2.1.2. Additionally, tramnet() needs a model matrix including the predictors, whose regression coefficients ought to be penalized. For numerical reasons it is useful to provide a scaled model matrix, instead of the original data, such that every parameter is equally affected by the regularization. Lastly, tramnet() will need the tuning parameters α ∈ [0, 1] and λ ∈ R + , with α representing a mixing parameter and λ controlling the extent of regularization. Setting λ = 0 will result in an unpenalized model, regardless of the value of α.
The R Journal Vol. XX/YY, AAAA 20ZZ ISSN 2073-4859 S3 methods accompanying the "tramnet" class will be discussed in section 2.3.

Censoring and likelihood forms
Specific combinations of F and the form of censoring yield log-log-concave log-likelihoods. Under these circumstances tramnet is not yet able to solve the resulting optimization problem. Table 1 indicates which model class can be fitted under what type of censoring in the current version of tramnet.

Prostate cancer data analysis
The regularized normal linear and extensions to transformed normal regression models will be illustrated using the Prostate data set (Stamey et al., 1989), which was used by Zou and Hastie (2005) to highlight properties of the elastic net.
R> data("Prostate", package = "lasso2") R> Prostate$psa <-exp(Prostate$lpsa) R> Prostate[, 1:8] <-scale(Prostate[, 1:8]) The data set contains 97 observations and 9 covariates. In the original paper the authors chose the log-transformed prostate specific antigen concentration (lpsa) as the response and used the eight remaining predictors log cancer volume (lcavol), log prostate weight (lweight), age of the patient (age), log benign prostatic hyperplasia amount (lbph), seminal vesicle invasion (svi coded as 1 for yes, 0 for no), log capsular penetration (lcp), Gleason score (gleason) and percentage Gleason score 4 or 5 (pgg45) as covariates. Zou and Hastie (2005) imposed an assumption on the conditional distribution of the response by log-transforming and fitting a linear model. In the following it is shown that the impact of this assumption may be assessed by estimating the baseline transformation from the data, followed by a comparison with the log-transformation applied by Zou and Hastie (2005). The linear models in lpsa and log(psa) are compared to transformation models with basis expansions in both log(psa) and psa, while specifying conditional normality of the transformed response. Additionally, the models are compared to an alternative implementation of regularized normal linear regression in penalized.

Linear and Box-Cox type regression models
Five different models will be used to illustrate important facettes of transformation models, including parametrization and interpretation. The models are summarized in Table 2 and will be elaborated on throughout this section. The comparison is based on unpenalized models first. Later, the section highlights the penalized models together with hyperparameter tuning.
R> fm_Pr <-psa~lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45 R> fm_Pr1 <-update(fm_Pr,~0 + .) R> x <-model.matrix(fm_Pr1, data = Prostate) The normal linear regression model is implemented in tram's Lm() function. Lm()'s parametrization differs from the usual linear model, hence caution has to be taken when interpreting the resulting Table 2: Summary of the five models illustrated in section 2.2, including their name throughout the manuscript, the R code to fit them and the mathematical formulation of their conditional cumulative distribution function. For comparison mp is included as an ordinary linear model, which is equivalent to model mt in terms of log-likelihood, but differs in the parametrization of the transformation function h and thus yields scaled coefficient estimates (cf. Table 3). Model mtp is a linear model parametrized in terms of a Bernstein basis of maximum order 1. This will yield the same coefficient estimates as mt but a log-likelihood that is comparable to models mt1 and mt2, whose transformation functions are parametrized in terms of higher order Bernstein bases. The log_first argument specifies whether the basis expansion is calculated on the log-transformed or untransformed response.
To recover the usual parametrization tramnet::coef.Lm() offers the as.lm = TRUE argument.
R> cfx_tramnet <-coef(mt, as.lm = TRUE) The transformation function for the linear model is depicted in Figure 1 (pink line). Because a linear baseline transformation imposes restrictive assumptions on the response's conditional distribution, it is advantageous to replace the linear baseline transformation by a more flexible one. In the case of the Box-Cox type regression model the linear baseline transformation h(y) = ϑ 1 + ϑ 2 log y is replaced by the basis expansion h(y) = a Bs,7 (log y) ϑ.
R> ord <-7 # flexible baseline transformation R> m01 <-BoxCox(psa~1, data = Prostate, order = ord, + extrapolate = TRUE, log_first = TRUE) R> mt1 <-tramnet(m01, x = x, alpha = 0, lambda = 0) The Box-Cox type regression model is then estimated with the BoxCox() function, while specifying the appropriate maximum order of the Bernstein polynomial. Because, the more flexible transformation slightly deviates from being linear, the normal linear model yields a smaller log-likelihood (cf. Table  3). To make sure that this improvement is not due to the increased number of parameters and hence overfitting, the models predictive capacities could be compared via cross-validation.
These results hold for the pre-specified log transformation of the response and a basis expansion thereof. Instead of prespecifying the log-transformation, its 'logarithmic nature' can be estimated from the data. Afterwards one can compare the deviation from a log-linear baseline transformation graphically and by inspecting the predictive performance of the model in terms of the out-of-sample log-likelihood.
R> m02 <-BoxCox(psa~1, order = 11, data = Prostate, extrapolate = TRUE) R> mt2 <-tramnet(m02, x = x, lambda = 0, alpha = 0) Indeed the baseline transformation in Figure 1 is similar to the basis expansion in the log-transformed response upon visual inspection. Because mt is estimated using the log-transformed response and mt1 and mt2 are based on the original scale of the response, the resulting model log-likelihoods are not comparable. To overcome this issue, one can fit a Box-Cox type model with maximum order 1, as this results in a linear, but alternatively parametrized baseline transformation.
R> m0p <-BoxCox(psa~1, order = 1, data = Prostate, log_first = TRUE) R> mtp <-tramnet(m0p, x = x, lambda = 0, alpha = 0) Figure 1 plots the three distinct baseline transformations resulting from models mt, mt1 and mt2. The initial assumption to model the prostate specific antigen concentration linearily on the log-scale seems to be valid when comparing the three transformation functions. The linear transformation in lpsa used in mt and the basis expansion in log(psa) (mt1) are almost indistinguishable and yield very similar coefficient estimates, as well as log-likelihoods (cf. Table 3, mtp vs. mt1). The basis expansion in psa (mt2) is expected to be less stable due to the highly skewed untransformed response. This is reflected in Figure 1, where the baseline transformation deviates from being linear towards the bounds of the response's support. However, the log-linear behaviour of h was clearly captured by this model and further supports the initial assumption of conditional log-normality of the response. For the same reasons, the resulting log-likelihood of mt2 is smaller than for mt1 (Table 3). Taken together, this exemplary analysis highlights the flexibility and usefulness of transformation models to judge crucial modelling assumptions. The original analysis prespecified a logtransformation of the response and then assumed conditional normality on this scale. Hence the baseline transformation of mt is of the form: h(lpsa) = ϑ 1 + ϑ 2 · lpsa. Now one can allow a more flexible transformation function in log(psa) to judge any deviations of h(log(psa)) from linearity, leading to a baseline transformation in terms of basis functions: a Bs,7 (log(psa)) ϑ in mt1. Lastly, instead of presuming a log-transformation one could estimate the baseline transformation from the raw response (psa), i.e. h(psa) = a Bs,11 (psa) ϑ in mt2. In this case, a higher order basis expansion was chosen to account for the skewness of the raw response. Notably, all three baseline transformations are fairly comparable. The basis expansion in psa deviates from being log-linear towards the boundaries of the response's support, as there are only few observations.

Hyperparameter tuning
This section features cross-validation, model-based optimization and profiling functions for hyperparameter tuning, whose appropriate values are highly problem-dependent and hard to know in advance. tramnet implements naive grid search and model-based optimization in the functions cvl_tramnet() and tramnet_mbo(), respectively. In the framework of regularized transformation models it is very natural to choose the out-of-sample log-likelihood as the objective function, because the notion of a mean square loss does not make sense for survival, let alone censored outcomes. The out-of-sample log-likelihood is in fact the log score, which is a proper scoring rule (Gneiting and Raftery, 2007). Comparison of the three transformation models on the Prostate data. Coefficient estimates are shown for each model, together with the in-sample log-likelihood in the last column. The first three models, mp, mt and mtp use a linear baseline transformation in lpsa and log(psa), respectively. The mp model was fit using penalized and gives the scaled version of the regression coefficients in mt, but the same log-likelihood. At the same time, mt and mtp differ only in their response variable and its subsequent log-transformation in mtp, yielding the same coefficient estimates but a different log-likelihood. Models mt1 and mt2 allow a flexible basis expansion in log(psa) and psa, respectively. Model mt1, allowing for a flexible basis expansion in lpsa, fits the data the best, however the resulting coefficient estimates are similar for all models. tramnet offers cross-validation in cvl_tramnet(), comparable to the optL1() and optL2() functions in penalized, which takes a sequence of values for λ and α and performs a simple -and arguably slow -grid search. Per default it computes 2-fold cross-validation, the user is encouraged, however, to judge the resulting bias-variance trade-off accordingly.
R> pen_cvl <-optL1(response = lpsa, penalized = x, lambda2 = 0, data = Prostate, + fold = 2) R> cvlt <-cvl_tramnet(object = mt, lambda = lambdas, alpha = 1, + folds = pen_cvl$fold) The resulting object is of class "cvl_tramnet" and contains a table for the cross-validated log-likelihoods for each fold and the sum thereof, the 'optimal' tuning parameter constellation which resulted in the largest cross-validated log-likelihood, tables for the cross-validated regularization paths, the folds and lastly the full fit based on the 'optimal' tuning parameters. Additionally, the resulting object can be used to visualize the log-likelihood and coefficient trajectories. These trajectories highly depend on the chosen folds and the user is referred to the full profiling functions discussed in section 2.2.2.

Model-based optimization
In contrast to naive grid search, model-based optimization comprises more elegant methods for hyperparameter tuning. tramnet offers the mbo_tramnet() and mbo_recommended() functions. The former implements Kriging-based hyperparameter tuning for the elastic net, the LASSO and ridge regression. mbo_tramnet() takes a "tramnet" object as input and computes the cross-validated loglikelihood based on the provided fold or folds arugment. The initial design is a random latin hypercube design with n_design rows per parameter. The number of sequential fits of the surrogate models is specified through n_iter and the range of the tuning parameters can be specified by max/min arguments. The default infill criterion is expected improvement. However, Bischl et al. (2017) encourage the use of the lower confidence bound over expected improvement, which can be achieved in mbo_tramnet() by specifying opt_crit = makeMBOInfillCritCB(). 10-fold cross-validation is used to compute the objective function for the initial design and each iteration. The recommended model is then extracted using mbo_recommended().
R> tmbo <-mbo_tramnet(mt, obj_type = "elnet", fold = 10) R> mtmbo <-mbo_recommended(tmbo, m0, x) Unlike in the previous section, one can directly optimize the tuning parameters in an elastic net problem instead of optimizing over one hyperparameter at a time or having to specify LASSO or ridge regression a priori. The output of mbo_tramnet() is quite verbose and can be shortened by using the helper function print_mbo().

R> print_mbo(tmbo)
## Recommended parameters: ## lmb=1.04e-05; alp=0.751 ## Objective: y = 710 Interpreting the output, model-based optimization suggests an unpenalized model with α = 0.75 and λ = 0. This result stresses the advantages of model-based optimization over naive or random grid search in terms of complexity and computational efficiency. In the end, the proposed model is unpenalized and thus does not introduce sparsity in the regression coefficients.

Regularization paths
As discussed before, it may be useful to inspect the full regularization paths over the tuning parameters λ and α. Akin to the functions profL1() and profL2() in package penalized, tramnet offers prof_lambda() and prof_alpha(). Since these functions take a fitted model of class "tramnet" as input, which is internally updated it is crucial to correctly specify the other tuning parameter in the model fitting step. In the example to come, mt was fit using α = 1 and λ = 0 resulting in a LASSO penalty only when profiling over λ. The resulting profile is depicted in Figure 2.
R> pfl <-prof_lambda(mt) prof_lambda() takes min_lambda, max_lambda and nprof as arguments and internally generates an equi-spaced sequence from min_lambda to max_lambda on the log scale of length nprof. By default this sequence ranges from 0 to 15 and is of length 5.

Additional constraints
In some applications, the specification of additional constraints on the shift parameters β are of interest. Most commonly, either positivity or negativity for some or all regression coefficients is aimed at. In tramnet additional inequality constraints can be specified via the constraints argument, which are internally handled as constraints [[1]] %*% beta > constraints [[2]]. Hence, to specify the constraint of strictly positive regression coefficients, one would supply an identity matrix of dimension p for the left hand side and the zero p-vector for the left hand side, as done in the following example.

S3 Methods
Building on the S3 infrastructure of the packages mlt and tram, this package provides corresponding methods for the following generics: coef(), logLik(), plot(), predict(), simulate() and residuals(). The methods' additional "tramnet"-specific arguments will be briefly discussed in this section.
coef.tramnet() suppresses the baseline transformation's coefficient estimatesθ by default and considers shift parameter estimatesβ below 10 −6 as 0 to stress the selected variables only. This threshold can be controlled by the tol argument. Hence, coef(mt,with_baseline = TRUE,tol = 0) returns all coefficients. Lastly, residuals.tramnet() computes the generalized residual r defined as the score contribution for sample i with respect to a newly introduced intercept parameter γ, which is restricted to be zero. In particular, yields the generalized residual with respect to γ for the model R> residuals(mtmbo)[1:5] ## [1] -6.50 -6.36 -6.60 -6.57 -4.17 In residual analysis and boosting it is common practice to check for associations between residuals and covariates that are not included in the model. In the prostate cancer example one could investigate, whether the variables age and lcp should be included in the model. To illustrate this particular case a non-parametric independence_test() from package coin can be used (Hothorn et al., 2008). First, the uncoditional transformation model m0 is fit. Afterwards, the tramnet models excluding age and lcp are estimated and their residuals extracted using the residuals.tramnet() method. Lastly, an independence test using a maximum statistic (teststat = "max") and a Monte Carlo based approximation Note that by default the plot is produced for each row in the design matrix. In unstratified linear transformation models this leads to shifted versions of the same curve on the transformation function's scale. A: Estimated conditional distribution function for every observation. B: Estimated conditional survivor function for every observation. The conditional survivor function is defined as S(y|x) = 1 − F Y (y|x). C: Conditional most likely transformation for every observation. Note that every conditional transformation function is a shifted version of the same curve. D: The conditional density for every observation can be calculated using f Y (y|x) = F (a(y) ϑ − x β)a (y) ϑ. E: A distribution function is fully characterized by its hazard function λ(y|x) = f Y (y|x)/S(y|x), which is depicted in this panel. F: The newdata argument can be used to plot the predicted most likely transformation for the provided data, in this case the first row of the Prostate data. of the null distribution based on resampling 10 6 times (distribution = approximate(1e6)), yields the results printed below. R> library("coin") R> m0 <-BoxCox(lpsa~1, data = Prostate, extrapolate = TRUE) R> x_no_age_lcp <-x[, !colnames(x) %in% c("age", "lcp")] R> mt_no_age_lcp <-tramnet(m0, x_no_age_lcp, alpha = 0, lambda = 0) R> r <-residuals(mt_no_age_lcp) R> it <-independence_test(r~age + lcp, data = Prostate, + teststat = "max", distribution = approximate(1e6)) R> pvalue(it, "single-step") ## age 0.023748 ## lcp <0.000001 Because there is substantial evidence against independence of the models' residuals to either lcp or age, we can conclude that it is worthwhile to include age and lcp in the model. Packages trtf (Hothorn, 2019b) and tbm (Hothorn, 2020a(Hothorn, , 2019a) make use of this definition of a residual for estimating and boosting transformation models, trees and random forests. For more theoretical insight the reader is referred to the above mentioned publications.