Motivation of Equation @ref(eq:cml-likelihood-nw).

It is instructive to start by considering an untruncated, censored observation where l = −∞, u = ∞ and m < v. The only information we obtain from the observation (m, v, l, u) is then that Y ∈ (m, v]. For deriving the relevant likelihood contribution, we may follow the stochastic approach to interval censored observations described in : let (C₁, C₂) denote a random vector in ℝ² that is independent of the target variable Y and which satisfies P(C₁ < C₂) = 1. Let

and define new random variables (M, V) = f(Y, C₁, C₂) by

Note that D can be reconstructed from (M, V): we have D = 0 if M = −∞, D = 1 if −∞ < M < V < ∞ and D = 2 if V = ∞.

It is instructive to proceed with the case where (C₁, C₂) and hence (M, V) is discrete. Then, for (m, v) ∈ supp(M, V) ∩ ℝ², we have

Likewise, for (m, v) ∈ supp(M, V) ∩ ({−∞} × ℝ), we obtain

and finally, for (m, v) ∈ supp(M, V) ∩ (ℝ × {∞}),

If we assume that the distribution of the censoring variable (C₁, C₂) is non-informative, i.e., its distribution does not depend on θ, the likelihood of observing (M, V) = (m, v) is equal to F_θ((m, v]), up to a factor that does not depend on θ. A similar argumentation can be used in the non-discrete case. Overall, noting that F_∞((−∞, ∞]) = 1, we have motivated the likelihood contribution F_θ((m, v]) ⋅ 1(m < v) for a censored, untruncated observation in @ref(eq:cml-likelihood-nw).

Next, consider an uncensored, truncated observation (m, v, l, u) where y = m = v; we may hence identify such an observation with (y, l, u). We may then proceed as in and assume that (L, U) is independent of Y and satisfies L ≤ U, with L possibly equal to −∞ and U possibly equal to ∞. Further, (L, U) shall have a density f_(L, U) with respect to some dominating sigma-finite measure ν. Truncation means that we only happen to observe (Y, L, U) if L < Y ≤ U. As a consequence, any observed value with M = V can be regarded as being drawn from the (μ ⊗ ν)-density

Subsequently, we write (Y^(t), L^(t), U^(t)) for a random vector following the above density, i.e.,

Conditioning this density on (L^(t), U^(t)) = (l, u), we arrive at an expression that does not involve the nuisance density f_(L, U):

Overall, we arrive at the (conditional) log-likelihood contribution log f_θ(y) − log F_θ((l, u]) for an uncensored, truncated observation in @ref(eq:cml-likelihood-nw).

Finally, truncation and censoring can occur at the same time, i.e., we have l ≤ m < v ≤ u with either l ≠ −∞ or u ≠ ∞. In accordance with the previous two cases, we make the assumption that Y, (C₁, C₂) and (L, U) are mutually independent and satisfy C₁ < C₂ and L < U. Define and

For simplicity, assume that all random variables are discrete. For any observation (m, v, l, u), one of the following four cases is met

In case l < m < v < u, we have

by the independence assumption. The factor in front does not depend on θ and is irrelevant for the (conditional) likelihood contribution. Likewise, in case l = m < v < u, we have

By definition of (M, V), the event in the numerator is the disjoint union of the following two sets:

By independence, we obtain that

Again, the factor in front of the fraction is independent of θ and is irrelevant for the likelihood. The two cases l < m < v = u and l = m < v = u can be treated similarly; in all cases, the likelihood contribution is equal to F_θ((m, v])/F_θ((l, u]) times a factor that does not depend on θ.

Related packages

For the less general cases of non-informative censoring without random truncation and fixed truncation, i.e., (L, U) constant for all observations, as well as for estimation of distribution parameters in the absence of censoring or random truncation, there are a number of R packages that can fit distributions, some of them also supporting weights. Among these are MASS (Venables and Ripley 2002), fitdistrplus (Delignette-Muller and Dutang 2015), survival (Therneau 2023), flexsurv (Jackson 2016). Note that fixed truncation is an operation that can be baked into the distribution family whose parameters are to be estimated, allowing for classical maximum likelihood estimation. Many of the packages also support classic regression of expected values given predictors. Distributional regression packages, such as gamlss (Stasinopoulos and Rigby 2007) and deepregression (Rügamer et al. 2023) currently do not support interval censoring or random truncation. See the following table for an overview of available features for each package.

Another R6-based interface is provided by ROOPSD (Robin 2022).

reservr builds upon the R packages tensorflow (Allaire and Tang 2022) and keras (Chollet, Allaire, et al. 2017) as an interface to the machine learning library TensorFlow (Abadi et al. 2015) to perform distributional regression. This underlying infrastructure is shared with the distributional regression package deepregression (Rügamer et al. 2023). The latter also supports distributional regression, but at the time of writing requires complete samples and does not support truncation or censoring.

The remaining parts of this paper are structured as follows: Section @ref(pkg-overview) details the core functionality of the corresponding R package reservr. It is split into definition of samples ℑ (Section @ref(trunc-obs)), definition of distribution families (Section @ref(distributions)), mathematical definitions of some available distribution families (Section @ref(dist-definitions)), estimation of distribution parameters (Section @ref(fit-dist)) and distributional regression using tensorflow (Section @ref(tensorflow)). A conclusion is given in Section @ref(conclusion).

Working with samples

A sample ℑ = {(m, v, l, u, w)_i} is represented as a tibble (from package tibble). The core function to create this tibble is trunc_obs(). A tibble created by trunc_obs() consists of five columns:

• x: If observed, the exact value of the random variable, referred to as Y in Section @ref(introduction). Otherwise NA.
• xmin: Lower interval censoring bound (M in Section @ref(introduction)) for the observation. If the observation is not censored, xmin is equal to x.
• xmax: Upper interval censoring bound (V in Section @ref(introduction)) for the observation. If the observation is not censored, xmax is equal to x.
• tmin: Lower truncation bound (L in Section @ref(introduction)). Only observations with x ≥ tmin are observed. Can be −∞ to indicate no lower truncation.
• tmax: Upper truncation bound (U in Section @ref(introduction)). Only observations with x ≤ tmax are observed. Can be ∞ to indicate no upper truncation.
• w: The weight associated with the observation. Defaults to 1.

Note that, unlike in Section @ref(introduction), the lower bounds of intervals in trunc_obs are included, that is, we allow for x ≥ tmin rather than x > tmin, and that the unknown variable of interest is called x instead of Y. For continuous random variables, the formulas are equivalent to the half-open formulation. For discrete random variables, xmin and tmin may have to be appropriately shifted, e.g., by replacing xmin by xmin − 0.5 for integer valued variables. The following code defines a sample of size 1 without truncation and censoring, with the realized value of 1.3.

trunc_obs(1.3)

## # A tibble: 1 × 6
##       x  xmin  xmax  tmin  tmax     w
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1   1.3   1.3   1.3  -Inf   Inf     1

Simulating randomly truncated and interval censored data from a standard normal distribution with 80% of the observations randomly interval censored and random uniform truncation L ∼ Unif[−2, 0] and U ∼ Unif[0, 2] can be simulated as follows

set.seed(123)
N <- 1000L
x <- rnorm(N)
is_censored <- rbinom(N, size = 1L, prob = 0.8) == 1L

c_lower <- runif(sum(is_censored), min = -2.0, max = 0.0)
c_upper <- c_lower + runif(sum(is_censored), min = 0, max = 1.0)

x_lower <- x
x_upper <- x

x_lower[is_censored] <- dplyr::case_when(
  x[is_censored] <= c_lower ~ -Inf,
  x[is_censored] <= c_upper ~ c_lower,
  TRUE ~ c_upper
)
x_upper[is_censored] <- dplyr::case_when(
  x[is_censored] <= c_lower ~ c_lower,
  x[is_censored] <= c_upper ~ c_upper,
  TRUE ~ Inf
)

t_lower <- runif(N, min = -2.0, max = 0.0)
t_upper <- runif(N, min = 0.0, max = 2.0)

is_observed <- t_lower <= x & x <= t_upper

obs <- trunc_obs(
  xmin = pmax(x_lower, t_lower)[is_observed],
  xmax = pmin(x_upper, t_upper)[is_observed],
  tmin = t_lower[is_observed],
  tmax = t_upper[is_observed]
)

Observations look like:

obs[8L:12L, ]

## # A tibble: 5 × 6
##        x    xmin   xmax   tmin  tmax     w
##    <dbl>   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 NA     -0.479   1.15  -1.93  1.15      1
## 2 NA     -0.177   1.79  -0.210 1.79      1
## 3 -0.556 -0.556  -0.556 -0.957 0.791     1
## 4 NA     -0.379   0.616 -0.379 0.616     1
## 5 NA      0.0575  1.45  -0.437 1.45      1

The total number of observations is smaller than the base population of 1000 due to truncation:

nrow(obs)

## [1] 623

The total number of censored observations is roughly 0.8 ⋅ nrow(obs).

sum(is.na(obs$x))

## [1] 496

In addition to the trunc_obs() constructor function, there are functions as_trunc_obs() for coercion, truncate_obs() for artificially changing truncation bounds, and repdel_obs() for computing randomly truncated reporting delay observations from general insurance claims data containing accident date, reporting delay and evaluation date information. The latter takes inputs of the form (T_acc, D, τ) where T_acc < τ are accident dates with corresponding reporting delays D ≥ 0 and τ is the calendar date of observation. It returns the sample (xmin = xmax = D, tmin = 0, tmax = τ − T_acc, w = 1) suitable for estimating the reporting delay distribution where a claim is only observed if it has been reported by the evaluation date, i.e., T_acc + D ≤ τ. Such an analysis was performed using reservr in .

Definition of distribution families

Distribution families are implemented using the R6 class system (Chang 2021). They inherit from the class Distribution and feature a common interface to

• manage fixed and free parameters of the underlying familiy,
• use basic distribution functions for random number generation and computation of the density, cumulative distribution, hazard and quantile function,
• use additional functions supporting parameter estimation procedures such as computing support or presence of a point mass,
• compile performance enhanced functions to speed up basic functions for repeated evaluation,
• provide tensorflow-specific implementations to support (deep) distributional regression.

A Distribution object represents a distribution family ℱ supported on a subset of the real line and parameterized by a fixed finite-dimensional parameter space Θ. The family may be a singleton, in which case it is rather a distribution than a distribution family.

reservr provides a set of basic distribution families, optionally with some fixed parameters, as well as transformations of distribution families that take one or more underlying distribution families. At the time of writing, these are:

dist <- dist_normal(sd = 1.0)

dist$sample(1L)

## Error in (function (n, mean = 0, sd = 1) : invalid arguments

set.seed(10L)
dist$sample(1L, with_params = list(mean = 0.0))

## [1] 0.01874617

set.seed(10L)
dist$sample(1L, with_params = list(mean = 0.0, sd = 2.0))

## [1] 0.03749234

set.seed(10L)
dist$sample(3L, with_params = list(mean = 0.0:2.0, sd = 0.5))

## [1] 0.009373085 0.907873729 1.314334725

dist <- dist_normal(sd = 1.0)
mix <- dist_mixture(dists = list(dist_normal(), NULL))

dist$default_params

## $mean
## NULL
## 
## $sd
## [1] 1

mix$default_params

## $dists
## $dists[[1]]
## A NormalDistribution with 2 dof
## 
## $dists[[2]]
## NULL
## 
## 
## $probs
## $probs[[1]]
## NULL
## 
## $probs[[2]]
## NULL

str(dist$get_placeholders())

## List of 1
##  $ mean: NULL

str(mix$get_placeholders())

## List of 2
##  $ dists:List of 2
##   ..$ :List of 2
##   .. ..$ mean: NULL
##   .. ..$ sd  : NULL
##   ..$ : NULL
##  $ probs:List of 2
##   ..$ : NULL
##   ..$ : NULL

str(dist$param_bounds)

## List of 2
##  $ mean:Classes 'Interval', 'R6' (-Inf, Inf) 
##  $ sd  :Classes 'Interval', 'R6' (0, Inf)

str(mix$param_bounds)

## List of 2
##  $ dists:List of 1
##   ..$ : NULL
##  $ probs:List of 1
##   ..$ :Classes 'Interval', 'R6' [0, 1]

str(dist$get_param_bounds())

## List of 1
##  $ mean:Classes 'Interval', 'R6' (-Inf, Inf)

str(mix$get_param_bounds())

## List of 2
##  $ dists:List of 1
##   ..$ :List of 2
##   .. ..$ mean:Classes 'Interval', 'R6' (-Inf, Inf) 
##   .. ..$ sd  :Classes 'Interval', 'R6' (0, Inf) 
##  $ probs:List of 2
##   ..$ :Classes 'Interval', 'R6' [0, 1] 
##   ..$ :Classes 'Interval', 'R6' [0, 1]

str(dist$get_param_constraints())

##  NULL

str(mix$get_param_constraints())

## function (params)

dist$get_components()

## list()

mix$get_components()

## [[1]]
## A NormalDistribution with 2 dof
## 
## [[2]]
## NULL

dist <- dist_normal()
flatten_params_matrix(dist$get_placeholders())

##      mean sd
## [1,]   NA NA

denscmp <- dist$compile_density()

if (requireNamespace("bench", quietly = TRUE)) {
  bench::mark(
    dist$density(-2:2, with_params = list(mean = 0.0, sd = 1.0)),
    denscmp(-2:2, matrix(c(0.0, 1.0), nrow = 5L, ncol = 2L, byrow = TRUE)),
    dnorm(-2:2, mean = rep(0.0, 5L), sd = rep(1.0, 5L))
  )
}

## # A tibble: 3 × 6
##   expression                            min  median `itr/sec` mem_alloc `gc/sec`
##   <bch:expr>                        <bch:t> <bch:t>     <dbl> <bch:byt>    <dbl>
## 1 dist$density(-2:2, with_params =… 26.09µs 28.59µs    34366.        0B     17.2
## 2 denscmp(-2:2, matrix(c(0, 1), nr…  4.06µs  4.83µs   201552.        0B      0  
## 3 dnorm(-2:2, mean = rep(0, 5L), s…  1.67µs  1.96µs   485721.        0B     48.6

Special families

Some of the distribution families available in reservr have tailored algorithms for parameter estimation, or are not commonly known. This section contains mathematical definitions of those function families.

dist1 <- dist_normal(mean = -1.0, sd = 1.0)
dist2 <- dist_exponential(rate = 1.0)
distb <- dist_blended(
  dists = list(dist1, dist2),
  breaks = list(0.0),
  bandwidths = list(1.0),
  probs = list(0.5, 0.5)
)

distt1 <- dist_trunc(dist1, min = -Inf, max = 0.0)
distt2 <- dist_trunc(dist2, min = 0.0, max = Inf)

distb1 <- distb$clone()
distb1$default_params$probs <- list(1.0, 0.0)
distb2 <- distb$clone()
distb2$default_params$probs <- list(0.0, 1.0)

Methods of estimating distribution parameters

This section describes the functions for the problem of estimating a parameter θ ∈ Θ given a sample ℑ and a parameterized family ℱ = {F_θ ∣ θ ∈ Θ}. Sometimes, the conditional log-likelihood in @ref(eq:cml-likelihood) can be directly maximized, yielding an estimate for θ. This is the default behavior in reservr if no specialized estimation routine for the provided family ℱ_θ is defined. Depending on whether there are box constraints, nonlinear constraints or no constraints on the parameter space Θ, different implementations of nonlinear optimization algorithms from nloptr (Johnson 2007), in particular truncated Newton (Dembo and Steihaug 1983) for unconstrained families, L-BFGS (Liu and Nocedal 1989) for box-constrained families and SLSQP (Kraft 1994) for general constrained families are employed.

In addition to the naive direct optimization approach, some families lend themselves to specialized estimation algorithms which usually show faster convergence due to making use of special structures in the parameter space Θ.

Estimating distribution parameters from truncated observations is handled using the generic fit() method. It delegates to fit_dist(), which is also generic with signature:

• dist: The distribution family to be fit
• obs: The trunc_obs object, or a vector of observed values
• start: Starting parameters, as a list compatible with dist$get_placeholders().

At the time of writing there are specialized algorithms for six types of families:

1. Blended distribution families
2. Erlang mixture distribution families
3. Generalized pareto distribution families with free lower bound u (estimated by the minimum of xmin over the sample)
4. Mixture distribution families
5. Translated distribution families with fixed offset and multiplier (transform the sample via $\tfrac{\cdot-\text{offset}}{\text{multiplier}}$ and fit the component distribution family to the transformed sample)
6. Uniform distribution families with free lower bound min or upper bound max (estimated by the minimum of xmin, for min, and the maximum of xmax, for max, over the sample)

If not present, the start parameter is obtained via the fit_dist_start() generic. This generic implements a family specific method of generating valid starting values for all placeholder parameters. A notable implementation is fit_dist_start.ErlangMixtureDistribution() for Erlang mixture distribution families. If the shape parameters are free, there are different initialization strategies that can be chosen using additional arguments to fit_dist_start():

• init = "shapes" paired with shapes = c(...) manually specifies starting shape parameters α
• init = "fan" paired with spread = d uses α = (1, 1 + d, …, 1 + (k − 1) ⋅ d) with a default of d = 1 resulting in α = (1, …, k)
• init = "kmeans" uses 1-dimensional K-means based clustering of the sample observations such that each cluster corresponds to a unique shape
• init = "cmm" uses the centralized method of moments procedure described in

Re-using dist <- dist_normal(sd = 1.0) from above and the generated sample obs, we can fit the free parameter mean:

dist <- dist_normal(sd = 1.0)
the_fit <- fit(dist, obs)
str(the_fit)

## List of 3
##  $ params:List of 1
##   ..$ mean: num 0.0822
##  $ opt   :List of 5
##   ..$ par        : Named num 0.0822
##   .. ..- attr(*, "names")= chr "mean"
##   ..$ value      : num 341
##   ..$ iter       : int 7
##   ..$ convergence: int 1
##   ..$ message    : chr "NLOPT_SUCCESS: Generic success return value."
##  $ logLik:Class 'logLik' : -341 (df=1)

Using the function plot_distributions() we can also assess the quality of the fit.

plot_distributions(
  true = dist,
  fitted = dist,
  empirical = dist_empirical(0.5 * (obs$xmin + obs$xmax)),
  .x = seq(-5, 5, length.out = 201),
  plots = "density",
  with_params = list(
    true = list(mean = 0.0, sd = 1.0),
    fitted = the_fit$params
  )
)

Here, the density labelled empirical corresponds to a kernel density estimate with automatic bandwidth selection.

We follow with an example of fitting an ErlangMixture(3) distribution family using various initialization strategies. Note that both, "kmeans" and "cmm" use the random number generator for internal K-means clustering. This necessitates setting a constant seed before running fit_dist_start() and fit() to ensure the chosen starting parameters are the same for both calls.

dist <- dist_erlangmix(list(NULL, NULL, NULL))
params <- list(
  shapes = list(1L, 4L, 12L),
  scale = 2.0,
  probs = list(0.5, 0.3, 0.2)
)

set.seed(1234)
x <- dist$sample(100L, with_params = params)

set.seed(32)
init_true <- fit_dist_start(dist, x, init = "shapes",
                            shapes = as.numeric(params$shapes))
init_fan <- fit_dist_start(dist, x, init = "fan", spread = 3L)
init_kmeans <- fit_dist_start(dist, x, init = "kmeans")
init_cmm <- fit_dist_start(dist, x, init = "cmm")
rbind(
  flatten_params(init_true),
  flatten_params(init_fan),
  flatten_params(init_kmeans),
  flatten_params(init_cmm)
)

##      shapes[1] shapes[2] shapes[3]    scale probs[1] probs[2] probs[3]
## [1,]         1         4        12 1.590800     0.43     0.33     0.24
## [2,]         1         4         7 2.688103     0.55     0.32     0.13
## [3,]         1         5        13 1.484960     0.43     0.36     0.21
## [4,]         2        10        24 1.010531     0.56     0.27     0.17

set.seed(32)
str(fit(dist, x, init = "shapes", shapes = as.numeric(params$shapes)))

## List of 4
##  $ params     :List of 3
##   ..$ probs :List of 3
##   .. ..$ : num 0.43
##   .. ..$ : num 0.33
##   .. ..$ : num 0.24
##   ..$ shapes:List of 3
##   .. ..$ : num 1
##   .. ..$ : num 4
##   .. ..$ : num 13
##   ..$ scale : num 1.59
##  $ params_hist: list()
##  $ iter       : int 1
##  $ logLik     :Class 'logLik' : -290 (df=6)

fit(dist, x, init = "fan", spread = 3L)$logLik

## 'log Lik.' -292.0026 (df=6)

fit(dist, x, init = "kmeans")$logLik

## 'log Lik.' -289.2834 (df=6)

fit(dist, x, init = "cmm")$logLik

## 'log Lik.' -293.1273 (df=6)

It should be noted that the different initialization methods had a considerable impact on the outcome in the example due to the discrete nature of Erlang mixture distribution shape parameters and thus the combinatorial difficulty of picking optimal shapes α. The fit() result for Erlang mixture distribution families contains an element named "params_hist". This can be populated by passing trace = TRUE to fit() and will record parameters after all ECME steps in the ECME-based estimation algorithms from and . The element "iter" contains the number of full ECME-Iterations that were performed.

Distributional regression using tensorflow integration

The maximization problem @ref(eq:cml-regression-likelihood) is delegated to tensorflow, which supplies ample stochastic optimization algorithms. Functions in reservr are necessary to create a suitable output layer for tensorflow that maps onto Θ and to provide an implementation of the (negative) log-likelihood in @ref(eq:cml-regression-likelihood) as a loss function. These two tasks are combined in tf_compile_model(). The function returns an object of class reservr_keras_model, which can be used for the estimation procedure.

Given input layers inputs and an intermediate output layer intermediate_output as well as a family of distributions dist, the function

• Compiles the loss for dist defined by @ref(eq:cml-regression-likelihood) as $l(g) = -\tfrac{1}{\#(\mathfrak{I}_{\mathrm{reg}})} \ell(g|\mathfrak{I}_{\mathrm{reg}})$, optionally disabling censoring or truncation for efficiency.
• Creates a list of final output layers mapping intermediate_output onto the parameter space Θ of dist using Distribution$tf_compile_params(). This step adds additional degrees of freedom to the overall model, and the approach is described in
• Runs keras3::compile() on the underlying keras.src.models.model.Model.

The following example defines a linear model with homoskedasticity assumption and fits it using 100 iterations of the Adam optimization algorithm (Kingma and Ba 2015). First, we simulate data (Y, X) from the model defined by X ∼ Unif(10, 20) and Y|X = x ∼ 𝒩(μ = 2x, σ = 1).

set.seed(1431L)
keras3::set_random_seed(1432L)

dataset <- tibble::tibble(
  x = runif(100, min = 10, max = 20),
  y = 2 * x + rnorm(100)
)

Next, we specify the distribution family ℱ = {𝒩(μ, σ = 1)|μ ∈ ℝ}, incorporating the homoskedasticity assumption.

dist <- dist_normal(sd = 1.0)

Using keras, we define an empty neural network, just taking x as an input and performing no transformation.

nnet_input <- keras3::keras_input(shape = 1L, name = "x_input")
nnet_output <- nnet_input

Then, tf_compile_model() adapts the input layer to the free parameter space Θ = ℝ. This introduces two parameters to the function family 𝒢 and implies the functional relationship μ = g(x) := θ₁ ⋅ x + θ₀. Since our sample is fully observed, we disable censoring and truncation, leading to the simplified loss

where f_μ(y) is the density of 𝒩(μ = μ, σ = 1) evaluated at y.

nnet <- tf_compile_model(
  inputs = list(nnet_input),
  intermediate_output = nnet_output,
  dist = dist,
  optimizer = keras3::optimizer_adam(learning_rate = 0.1),
  censoring = FALSE,
  truncation = FALSE
)
nnet$dist
nnet$model

The fit can now be performed, modifying the parameters (weights) of nnet in-place. Note that the argument y of fit accepts a trunc_obs object. In our example, the vector y is silently converted to an untruncated, uncensored trunc_obs object. fit() returns the keras_training_history of the underlying call to fit() on the keras.src.models.model.Model.

nnet_fit <- fit(
  nnet,
  x = dataset$x,
  y = dataset$y,
  epochs = 100L,
  batch_size = 100L,
  shuffle = FALSE,
  verbose = FALSE
)

The training history can be plotted, displaying the loss by epoch (black circles), and a blue smoothing line.

# Fix weird behavior of keras3
nnet_fit$params$epochs <- max(nnet_fit$params$epochs, length(nnet_fit$metrics$loss))
plot(nnet_fit)

The predict() method of reservr_keras_model takes input tensors and returns the predicted distribution parameters as a list compatible with dist$get_placeholders(). We can thus extract the only parameter mean and compare it to an OLS fit for the same dataset:

pred_params <- predict(nnet, data = list(keras3::as_tensor(dataset$x, keras3::config_floatx())))

lm_fit <- lm(y ~ x, data = dataset)

dataset$y_pred <- pred_params$mean
dataset$y_lm <- predict(lm_fit, newdata = dataset, type = "response")

library(ggplot2)
ggplot(dataset, aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = y_pred), color = "blue") +
  geom_line(aes(y = y_lm), linetype = 2L, color = "green")

Since a reservr_keras_model includes the underlying keras.src.models.model.Model, its parameters can also be extracted and compared to the OLS coefficients

coef_nnet <- rev(as.numeric(nnet$model$get_weights()))
coef_lm <- unname(coef(lm_fit))

str(coef_nnet)
str(coef_lm)

We now discuss a more complex example involving censoring, using the right-censored ovarian dataset bundled with the survival package (R Core Team 2023). Our goal is to predict the rate parameter of an exponential survival time distribution in cancer patients given four features X = (age, resid.ds, rx, ecog.ps) collected in the study. The variables resid.ds, rx and ecog.ps are indicator variables coded in {1, 2}. age is a continuous variable with values in (38, 75). Due to the different scale of the age variable, it is useful to separate it from the other variables in order to perform normalization. Normalization using keras3::layer_normalization() transforms its input variables to zero mean and unit variance. This step is not necessary for the categorical features.

set.seed(1219L)
tensorflow::set_random_seed(1219L)
keras3::config_set_floatx("float32")

dist <- dist_exponential()
ovarian <- survival::ovarian
dat <- list(
  y = trunc_obs(
    xmin = ovarian$futime,
    xmax = ifelse(ovarian$fustat == 1, ovarian$futime, Inf)
  ),
  x = list(
    age = keras3::as_tensor(ovarian$age, keras3::config_floatx(), shape = nrow(ovarian)),
    flags = k_matrix(ovarian[, c("resid.ds", "rx", "ecog.ps")] - 1.0)
  )
)

Next, we define the input layers and shapes, conforming to our input predictor list dat$x.

nnet_inputs <- list(
  keras3::keras_input(shape = 1L, name = "age"),
  keras3::keras_input(shape = 3L, name = "flags")
)

age will be normalized and then concatenated to the other features, stored in flags, resulting in a 4-dimensional representation. We then add two hidden ReLU-layers each with 5 neurons to the network and compile the result, adapting the 5-dimensional hidden output to the parameter space Θ = (0, ∞) for the rate parameter of an exponential distribution. This is accomplished using a dense layer with 1 neuron and the softplus activation function.

hidden1 <- keras3::layer_concatenate(
  keras3::layer_normalization(nnet_inputs[[1L]]),
  nnet_inputs[[2L]]
)
hidden2 <- keras3::layer_dense(
  hidden1,
  units = 5L,
  activation = keras3::activation_relu
)
nnet_output <- keras3::layer_dense(
  hidden2,
  units = 5L,
  activation = keras3::activation_relu
)

nnet <- tf_compile_model(
  inputs = nnet_inputs,
  intermediate_output = nnet_output,
  dist = dist,
  optimizer = keras3::optimizer_adam(learning_rate = 0.01),
  censoring = TRUE,
  truncation = FALSE
)
nnet$model

For stability reasons, the default weight initialization is not optimal. To circumvent this, we estimate a global exponential distribution fit on the observations and initialize the final layer weights such that the global fit is the initial prediction of the network.

str(predict(nnet, dat$x))
global_fit <- fit(dist, dat$y)
tf_initialise_model(nnet, params = global_fit$params, mode = "zero")
str(predict(nnet, dat$x))

Finally, we can train the network and visualize the predictions.

nnet_fit <- fit(
  nnet,
  x = dat$x,
  y = dat$y,
  epochs = 100L,
  batch_size = nrow(dat$y),
  shuffle = FALSE,
  verbose = FALSE
)

nnet_fit$params$epochs <- max(nnet_fit$params$epochs, length(nnet_fit$metrics$loss))
plot(nnet_fit)

ovarian$expected_lifetime <- 1.0 / predict(nnet, dat$x)$rate

A plot of expected lifetime by (age, rx) shows that the network learned longer expected lifetimes for lower age and for treatment group (rx) 2. The global fit is included as a dashed blue line.

Individual predictions and observations can also be plotted on a subject level.