---
abstract: 'We introduce a general class of continuous univariate distributions with positive support obtained by transforming the class of two-piece distributions. We show that this class of distributions is very flexible, easy to implement, and contains members that can capture different tail behaviours and shapes, producing also a variety of hazard functions. The proposed distributions represent a flexible alternative to the classical choices such as the log-normal, Gamma, and Weibull distributions. We investigate empirically the inferential properties of the proposed models through an extensive simulation study. We present some applications using real data in the contexts of time-to-event and accelerated failure time models. In the second kind of applications, we explore the use of these models in the estimation of the distribution of the individual remaining life.'
author:
- '[Francisco J. Rubio[^1]]{} and Yili Hong[^2]'
title: Survival and lifetime data analysis with a flexible class of distributions
---

[*Key Words: AFT model; composite models; hazard function; logarithmic transformation; remaining life..*]{}

Introduction {#sec:intro}
============

In many areas, including medical applications, the quantities of interest take positive values. For instance, in survival analysis, the interest typically focuses on modelling the survival times of a group of patients in terms of a set of covariates (see *e.g.* [@L03]). Other areas where positive observations appear naturally are finance (*e.g.* in modelling the size of reinsurance claims), network traffic modelling [@M01], reliability theory [@ME98], environmental science [@MG10], among many others. Parametric distributions provide a parsimonious way of describing the distribution of those quantities. Some of the most popular choices for modelling positive observations are the lognormal, log-logistic, Gamma, and Weibull distributions. We refer the readers to [@M12] for an extensive overview of these sorts of distributions as well as a study of their inferential properties in the presence of censored observations. However, these distributions do not always provide a good fit of the data. For example, when the data present heavier tails and/or a different shape around the mode than those captured by these distributions. In recent years, there has been an increasing interest in the development of flexible distributions with positive support in order to cover departures from the classical choices. Two popular strategies for generating new flexible distributions with positive support consist of:

*Adding a shape parameter to an existing distribution with positive support*. For instance, in the context of reliability and survival analysis, [@MO97] proposed a transformation of a distribution $F(y;\theta)$, $y>0$, that introduces a new parameter $\gamma>0$. This transformation is defined through the cumulative distribution function (CDF) $$\begin{aligned}
\label{MOdistribution}
G(y;\theta,\gamma)=\dfrac{ F(y;\theta)}{F(y;\theta)+\gamma(1-F(y;\theta))}.\end{aligned}$$

The interpretation of the parameter $\gamma$ is given in [@MO97] in terms of the behavior of the ratio of hazard rates of $F$ and $G$. This ratio is increasing in $y$ for $\gamma\geq 1$ and decreasing in $y$ for $0<\gamma\leq 1$. This transformation is then proposed for the exponential and Weibull distributions in [@MO97] in order to generate more flexible models for lifetime data. Clearly, for $\gamma=1$, $G$ and $F$ coincide. Many choices of $F(;\theta)$ have already been studied in the literature. We refer the reader to [@FS06] for a general mechanism for adding parameters to a distribution.

*Using transformations from ${\mathbb R}$ to ${\mathbb R}_+$*. The most common choice for this transformation is the exponential function. The idea is to define a positive variable $Y$ by transforming a real variable $X$ through $Y=\exp(X)$. This method is used to produce the class of log-symmetric distributions. This is, the family of positive random variables such that their logarithm is symmetrically distributed. Some members of this class are the lognormal, log-logistic, and log-Student-$t$ distributions, which are obtained by transforming the normal, logistic, and Student-$t$ distributions, respectively. More recently, other families of distributions have been proposed by using this idea, such as the log Birnbaum-Saunders distribution [@B08], log skew-elliptical distributions [@MG10], log-generalised extreme value distributions [@RD15], and log-scale mixtures of normals [@VS15].

In this paper, we propose a new class of flexible distributions with support on ${\mathbb R}_+$ by applying the second method to the family of two-piece distributions [@FS98; @A05; @RS14]. In Section \[LTPDistributions\], we introduce the proposed class of distributions and show that it contains very flexible members that can capture a wide variety of shapes and tail behaviours. We show that these models can be seen as a subclass of composite models, which are of great interest in finance. The associated hazard functions are non-monotone with either increasing or decreasing right tails. These distributions are easy to implement using the R packages [@R13] ‘twopiece’ and ‘TPSAS’, which are available upon request. In Section \[MaximumLikelihoodEstimation\], we discuss the properties of the maximum likelihood estimators (MLEs) associated to these models. Although a formal study of the asymptotic properties of the proposed models is beyond the scope of this paper, we present a simulation study which reveals that adding a shape parameter, via two-piece transformations, has little effect on the performance of the maximum likelihood estimators. In Section \[NumericalExamples\], we present two kinds of applications with real data. In the first example, we illustrate the use of the proposed distributions in the context of data fitting. The main application is presented in the second and third examples, where we employ the proposed distributions for modelling the errors in an accelerated failure time model (AFT) with applications to medical data. In the third example, we discuss the use of a certain class of prediction intervals of the remaining life, which are informative for individual prognosis. In all of these examples, we discuss model selection between some appropriate competitors and the selection of the baseline distribution in the proposed family of distributions.

Log Two-piece Distributions {#LTPDistributions}
===========================

For the sake of completeness, let us first recall the definition of two-piece distributions. Let $s(\cdot;\delta)$ be a symmetric unimodal density, with mode at $0$, with support on $\mathbb R$, and let $\delta\in\Delta\subset{\mathbb R}$ be a shape parameter (location and scale parameters can be added in the usual way). The corresponding CDF will be denoted as $S(\cdot;\delta)$. The shape parameter $\delta$ typically controls the tails of the density. For example, in the cases where $s(\cdot;\delta)$ is either a Student-$t$ density with $\delta>0$ degrees of freedom or an exponential power density with power parameter $\delta>0$ (see Appendix A).

A real random variable $X$ is said to be distributed according to a two-piece distribution if its probability density function (PDF) is given by (see *e.g.* [@RS14]): $$\begin{aligned}
\label{TP}
s_{tp}(x;\mu,\sigma_1,\sigma_2,\delta) = \dfrac{2}{\sigma_1+\sigma_2}\left[s\left(\dfrac{x-\mu}{\sigma_1};\delta\right) I(x<\mu) + s\left(\dfrac{x-\mu}{\sigma_2};\delta\right) I(x\geq \mu) \right].\end{aligned}$$

This is, a two-piece density is obtained by continuously joining two half-$s$ densities with different scale parameters on either side of the location $\mu$. The density (\[TP\]) is unimodal, with mode at $\mu$, it is asymmetric for $\sigma_1\neq \sigma_2$, and coincides with the original density $s$ for $\sigma_1=\sigma_2$. Moreover, the tail behaviour of the PDF in (\[TP\]) is the same in each direction, by construction. A popular reparameterisation is obtained by redefining $\sigma_1=\sigma a(\gamma)$ and $\sigma_2 = \sigma b(\gamma)$, where $a(\cdot)$ and $b(\cdot)$ are positive functions of the parameter $\gamma$ [@A05]. Two common choices for $a(\cdot)$ and $b(\cdot)$ are the inverse scale factors $\{a(\gamma),b(\gamma)\}=\{\gamma,1/\gamma\}$, $\gamma \in {\mathbb R}_+$ [@FS98], and the epsilon-skew parameterisation $\{a(\gamma),b(\gamma)\}=\{1-\gamma,1+\gamma\}$, $\gamma \in (-1,1)$ [@MH00]. Other parameterisations are explored in [@RS14]. The PDF associated to this reparameterisation is given by $$\begin{aligned}
\label{TP2}
s_{tp}(x;\mu,\sigma,\gamma,\delta) = \dfrac{2}{\sigma[a(\gamma)+b(\gamma)]}\left[s\left(\dfrac{x-\mu}{\sigma b(\gamma)};\delta\right) I(x<\mu) + s\left(\dfrac{x-\mu}{\sigma a(\gamma)};\delta\right) I(x\geq \mu) \right]\end{aligned}$$ This transformation preserves the existence of moments and the ease of use of the original distribution $s$. The corresponding cumulative distribution function and quantile function can be easily obtained from this expression (see [@A05]). This class of distributions has been shown to have good inferential properties for regular choices of the baseline density $s$ [@A05; @JA10].

By applying method (ii), described in Section \[sec:intro\], to the family of two-piece distributions, we can produce distributions with support on ${\mathbb R}_+$ as follows.

A positive random variable $Y$ is said to be distributed according to a log two-piece (LTP) distribution if its PDF is given by: $$\begin{aligned}
\label{LTPPDF}
s_l(y; \mu,\sigma,\gamma,\delta) = \dfrac{2}{y\sigma[a(\gamma)+b(\gamma)]} &\left[ s\left(\dfrac{\log(y)-\mu}{\sigma b(\gamma)};\delta\right)I(y<e^{\mu})\right.\nonumber\\ &\quad\quad+\left.s\left(\dfrac{\log(y)-\mu}{\sigma a(\gamma)};\delta\right)I(y\geq e^{\mu}) \right].\end{aligned}$$

Given that the class of two-piece distributions contains all the symmetric unimodal distributions with support on the real line, it follows that the class of LTP distributions contains the class of log-symmetric distributions as well as models such that the distribution of $\log Y$ is asymmetric. The LTP Laplace distribution, which is obtained by using a Laplace baseline density $s$ in (\[LTPPDF\]), has been studied in [@K01]. However, other types of log two–piece distributions have not been studied to the best of our knowledge. The corresponding CDF is given by $$\begin{aligned}
\label{LTPCDF}
S_l(y; \mu,\sigma,\gamma,\delta) &=& \dfrac{2b(\gamma)}{a(\gamma)+b(\gamma)}  S\left(\dfrac{\log(y)-\mu}{\sigma b(\gamma)};\delta\right)I(y<e^{\mu}) \notag\\
&+&\left[\dfrac{b(\gamma)-a(\gamma)}{a(\gamma)+b(\gamma)} + \dfrac{2a(\gamma)}{a(\gamma)+b(\gamma)}S\left(\dfrac{\log(y)-\mu}{\sigma a(\gamma)};\delta\right)\right]I(y\geq e^{\mu}) .\end{aligned}$$

We can observe that the ratio of the mass cumulated on either side of the value $y=e^{\mu}$ is given by $$\begin{aligned}
R(\gamma) &=& \dfrac{S_l(e^{\mu}; \mu,\sigma,\gamma,\delta)}{1-S_l(e^{\mu}; \mu,\sigma,\gamma,\delta)} = \dfrac{b(\gamma)}{a(\gamma)}.\end{aligned}$$

This helps us to identify the different roles of the parameters $\gamma$ and $\delta$. The parameter $\gamma$ controls the cumulation of mass on either side of $y=e^{\mu}$, while the parameter $\delta$ controls the tails of the density. In Figure \[fig:LTPN\]a we present some examples of a two-piece normal PDF with different values of the parameter $\gamma$. In these cases, the parameter $\gamma$ only affects the asymmetry of the density. Figure \[fig:LTPN\]b shows the corresponding LTP normal PDFs. We can observe that in these cases the parameter $\gamma$ affects the shapes of the density. That is, it controls the mass cumulated above and below the value $y=1$ as well as the spread and mode of the density. The corresponding hazard function can be easily constructed from (\[LTPPDF\]) and (\[LTPCDF\]). Figure \[fig:DH\] shows the variety of shapes of the density and hazard functions obtained for a log two-piece sinh-arcsinh distribution (LTP SAS, which is obtained by using a symmetric sinh-arcsinh baseline density function in (\[TP2\]), see also [@R15]. The corresponding expression is provided in Appendix A). The implementation of LTP distributions is straightforward in R by using the packages ‘twopiece’ and ‘TPSAS’, which are freely available upon request. Moreover, the $p$th moment of a LTP distribution exists, whenever the $p$th moment of the underlying (log-symmetric) log-$s$ distribution exists. In particular, all moments of the LTP normal distribution exist.

  ----- -----
        
   (a)   (b)
  ----- -----

  ----- -----
        
   (a)   (b)
        
   (c)   (d)
        
   (e)   (f)
  ----- -----

An alternative construction {#an-alternative-construction .unnumbered}
---------------------------

The family of two-piece distributions (\[TP\])-(\[TP2\]) can be seen as a special kind of finite mixtures of truncated PDFs, as shown in [@RS14]. In a similar fashion, the family of log two-piece distributions can be obtained as a particular class of finite mixtures of truncated distributions with positive support. In the context of survival and size distributions these sorts of mixtures are known as composite models (see [@NB13] for a literature review). Recall first that the PDF of a composite model can be written as: $$\begin{aligned}
\label{CompPDF}
s_c(y) = \omega s_1^*(y) I(y\leq \theta) + (1-\omega) s_2^*(y) I(y>\theta),\,\,\, y>0,\end{aligned}$$ where $\omega = \dfrac{s_2(\theta)S_1(\theta)}{s_2(\theta)S_1(\theta) + s_1(\theta)[1-S_2(\theta)]}$, $\theta>0$ is a threshold parameter, $s_1^*(y)=\dfrac{s_1(y)}{S_1(\theta)}$, $s_2^*(y)=\dfrac{s_2(y)}{1-S_2(\theta)}$, $s_1$ and $s_2$ are continuous PDFs with support on ${\mathbb R}_+$, and $S_1$ and $S_2$ are the corresponding CDFs. If we fix $s_1(y)=\dfrac{1}{\sigma_1 y}s\left(\dfrac{\log(y)-\mu}{\sigma_1}\right)$, $s_2(y)=\dfrac{1}{\sigma_1 y}s\left(\dfrac{\log(y)-\mu}{\sigma_2}\right)$, for some symmetric density $s$ with support on ${\mathbb R}$, and $\theta=e^{\mu}$, then it follows that (\[CompPDF\]) coincides with (\[LTPPDF\]), up to a reparameterisation.

From this alternative construction, we conclude that the family of LTP distributions represents a subclass of composite models with the appealing properties and interpretability of parameters discussed above. This also allows us to motivate the use of LTP distributions as survival and size distributions.

Models and Maximum Likelihood Estimation {#MaximumLikelihoodEstimation}
========================================

In this section, we present the parameter estimation procedure for time-to-event and accelerated failure time (AFT) models.

Time-to-event model
-------------------

Let ${\bf T}=(T_1,\dots,T_n)$ be an independent sample of survival times distributed as in (\[LTPPDF\]). The likelihood function of the parameters $(\mu,\sigma,\gamma,\delta)$ is defined as: $$\begin{aligned}
L(\mu,\sigma,\gamma,\delta) = \prod_{j=1}^n s_l(T_j;\mu,\sigma,\gamma,\delta).\end{aligned}$$ The MLE is defined as the parameter values that maximise the likelihood function. By noting that $$\begin{aligned}
L(\mu,\sigma,\gamma,\delta) \propto \prod_{j=1}^n s_{tp}(\log(T_j);\mu,\sigma,\gamma,\delta),\end{aligned}$$ it follows that the MLEs of the parameters of LTP distributions are the same as the MLEs of the parameters of the underlying two-piece distribution for the sample $\log({\bf T}) = [\log(T_1),\dots,$ $\log(T_n)]$. Inferential aspects of 3- and 4-parameter two-piece distributions have been largely discussed. For example, [@A05] show that, under certain regularity conditions on the baseline density $s$ in (\[TP2\]), the maximum likelihood estimators of the parameters of these distributions are consistent and asymptotically normal under the epsilon-skew parameterisation. [@JA10] and [@RS14] study some parameterisations that induce parameter orthogonality between the parameters $\mu$ and $\sigma$, showing that the epsilon-skew parameterisation induces this property. Parameter orthogonality, in turn implies a good asymptotic behaviour of the MLE [@JA10]. In most cases, the MLE is not available in closed-form, and it has to be obtained numerically.

Samples containing censored observations are common in the context of survival analysis. The most common types of censoring in this context correspond to:

(i) Left-censoring: when the phenomenon of interest has already occurred before the start of the study. A left-censored observation is an interval of the type $[0,T_j)$, where $T_j$ represents the start of the study for subject $j$.

(ii) Interval censoring: when the phenomenon of interest occurs within a finite period of time $[T_j^L,T_j^R]$.

(iii) Right-censoring: when the phenomenon of interest is not observed during the duration of the study. A right-censored observation is an interval of the type $(T_j,\infty]$, where $T_j$ represents the duration of the study for subject $j$.

Ignoring censoring induces bias in the estimation of the parameters. Different types of censoring imply different contributions of the observations to the likelihood function. The contribution of a left-censored observation to the likelihood is $S_l(T_j;\mu,\sigma,\gamma,\delta)$; while the contribution of an interval-censored observation to the likelihood is $S_l(T_j^R;\mu,\sigma,\gamma,\delta)-S_l(T_j^L;\mu,\sigma,\gamma,\delta)$; and the contribution of a right-censored observation to the likelihood is $1-S_l(T_j;\mu,\sigma,\gamma,\delta)$. If we define the sets $\text{Left}=\{j: T_j \text{ is left-censored}\}$, $\text{Int}=\{j: T_j \text{ is interval-censored}\}$, $\text{Right}=\{j: T_j \text{ is right-censored}\}$, and $\text{Obs}=\{j: T_j \text{ in uncensored}\}$, then we can write the likelihood function as follows: $$\begin{aligned}
L(\mu,\sigma,\gamma,\delta) = &\prod_{j\in\text{Obs}} s_l(T_j;\mu,\sigma,\gamma,\delta) \times \prod_{j\in\text{Left}} S_l(T_j;\mu,\sigma,\gamma,\delta) \\&\times  \prod_{j\in\text{Right}}\left[1-S_l(T_j;\mu,\sigma,\gamma,\delta)\right]
\times  \prod_{j\in\text{Int}} \left[ S_l(T_j^R;\mu,\sigma,\gamma,\delta)-S_l(T_j^L;\mu,\sigma,\gamma,\delta)\right].\end{aligned}$$

The latter expression emphasises the practical importance of using distributions with a tractable distribution function.

Accelerated failure time models
-------------------------------

AFT models are a useful tool for modelling the set of survival times ${\bf T}=(T_1,\dots,T_n)$ in terms of a set of covariates $\bm{\beta} = (\beta_1,\dots,\beta_p)$ through the model equation: $$\begin{aligned}
\label{AFTModel}
h(T_j) = x_j^{\top}\bm{\beta} + \varepsilon_j, \,\,\, j=1,\dots,n,\end{aligned}$$

where ${\bf X}= (x_1,\dots,x_n)^{\top}$ is an $n\times p$ known design matrix and $\varepsilon_j \stackrel{ind.}{\sim} F(\cdot;\bm{\theta})$, $F$ is a continuous distribution with support on ${\mathbb R}$ and parameters $\bm{\theta}\in \Theta \subset {\mathbb R}^d$, and $h:{\mathbb R}_+\rightarrow{\mathbb R}$ is a continuous increasing function. The most common choice for $h$ is the logarithmic function, while the distribution of the errors $\varepsilon_j$ is typically assumed to be normal. Given that the assumption of normality of the errors can be restrictive in practice, other distributional assumptions have been recently studied such as the log Birnbaum-Saunders [@B08], finite mixtures of normal distributions [@KL08], the symmetric family of scale mixtures of normals [@VS15], and the log-generalised extreme value distribution [@RD15].

AFT models are extremely relevant in medicine, given that survival data naturally arise in many medical studies, which typically involve the follow-up of other covariates. The presence of different types of censored observations is common in this context [@B08; @KL08; @VS15]. If we assume that the errors $\varepsilon_j$ are distributed according to a LTP distribution with location $0$ and $\bm{\theta}=(\sigma,\gamma,\delta)$, then we can write the likelihood function as follows, $$\begin{aligned}
L(\bm{\beta},\sigma,\gamma,\delta) =& \prod_{j\in\text{Obs}} s_l(T_j;x_j^{\top}\bm{\beta},\sigma,\gamma,\delta) \times \prod_{j\in\text{Left}} S_l(T_j;x_j^{\top}\bm{\beta},\sigma,\gamma,\delta)\\ &\times  \prod_{j\in\text{Right}}\left[1-S_l(T_j;x_j^{\top}\bm{\beta},\sigma,\gamma,\delta)\right]\\
&\times  \prod_{j\in\text{Int}} \left[ S_l(T_j^R;x_j^{\top}\bm{\beta},\sigma,\gamma,\delta)-S_l(T_j^L;x_j^{\top}\bm{\beta},\sigma,\gamma,\delta)\right],\end{aligned}$$

with the notation discussed previously. It is important to notice that by using asymmetric errors, we obtain a curve that does not represent the mean response. However, as discussed in [@AG08], this lack of centring can be calibrated after estimating the parameters by adding a suitable quantity $M_{\varepsilon}$ which reflects the lack of centring of the errors. For instance, in order to obtain the mean response, we can use $M_{\varepsilon}=-{\mathbb E}[\varepsilon_j]$, computed at the MLE of the parameters of the error distribution. This strategy will only affect the intercept parameter. When using baseline models with infinite variance (such as a log-Cauchy distribution), one might opt for centring around the median (or another quantile), instead of the mean. A formal study of the asymptotic properties of the MLEs under different types of censoring is beyond the scope of this paper. However, in Section \[SimulationStudy\], we illustrate the performance of the MLEs in a linear regression model with censored observations through a simulation study.

Simulation Study {#SimulationStudy}
================

In this section, we present a simulation study in order to illustrate the performance of the MLEs of the parameters of some LTP distributions. Throughout, we employ the epsilon-skew parameterisation discussed previously. In our first simulation scenario, we simulate $N=10,000$ samples of sizes $n=30,50,100,250,500,1000$ from a LTP normal (log two-piece normal) with different combinations of the parameter values: $\mu=0$, $\sigma=1$, and $\gamma=0,0.25,0.5,0.75$. Negative values of $\gamma$ would induce similar results, since they produce the corresponding reflected density about $e^{\mu}$, and are therefore omitted. For each of these samples, we calculate the corresponding MLEs, using the R command ‘optim’, and calculate the bias, variance, and root-mean-square error (RMSE) of these. In our second simulation scenario, we simulate $N=10,000$ samples from a LTP $t$ (log two-piece Student-$t$) with parameters: $\mu=0$, $\sigma=1$, $\gamma=0,0.25,0.5,0.75$, and $\delta=1$. The third simulation scenario is analogous to the second scenario, with $\delta=2$. In the fourth scenario, we simulate $N=10,000$ samples from a LTP SAS with parameters : $\mu=0$, $\sigma=1$, $\gamma=0,0.25,0.5,0.75$, and $\delta=0.75$. Tables \[table:LTPN\]–\[table:LTPSAS2\] present the results of these simulations.

In second class of simulations, we investigate the performance of the use of log two-piece errors in AFT models (\[AFTModel\]). For this purpose, we simulate from the linear regression model: $$\begin{aligned}
\log(y_j) = {\bf x}_j^{\top}\bm{\beta} + \varepsilon_j, \,\,\, j=1,\dots,n,\end{aligned}$$ with $n=100,250,500$, $\beta = (1,2,3)^{\top}$, and ${\bf x}_j = (1,x_{j1},x_{j2})^{\top}$. The second and third entries of the covariates ${\bf x}_j$ are simulated from a right-half-normal with scale parameter $1/3$. For the distribution of the errors $\varepsilon_j$ we consider the following cases: (i) a two-piece normal distribution with parameters $\mu=0$, $\sigma=0.25$, and $\gamma=0,0.25,0.5$; and (ii) a TP SAS distribution [@R15] with parameters $\mu=0$, $\sigma=0.25$, $\gamma=0,0.25,0.5$, and $\delta=0.75$. We truncate the observations $y_j$ that are greater than $17.5$. This censoring mechanism produces samples with 15%–35% censored observations. Tables \[table:LTPNR\]–\[table:LTPSASR\] present the results of these simulations.

The overall conclusions of this extensive simulation study are that the value of the shape parameter $\gamma$ does not seem to greatly affect the performance of the MLEs, while the use of models with a tail parameter $\delta$ have a clear effect on the performance of the MLEs. The performance of the MLEs of $\delta$ in LTP $t$ and LTP SAS models for small samples is different: the bias is smaller in the LTP SAS model. However, the estimation of $\sigma$ is more accurate in the LTP $t$ model. This is, perhaps, an unsurprising conclusion, given that it is well–known that it is difficult to learn about tail parameters with small samples and that tail parameters control the tail behaviour differently in different models. However, this analysis helps us to quantify the order of observations required for an accurate estimation. For LTP models with 4 parameters, such as the LTP $t$ and LTP SAS models, it is necessary to have at least $200$ observations in order to accurately estimate the tail parameters. In fact, the proposed flexible models are not recommended with small samples since, intuitively, these do not contain information about the features captured by the shape parameters $\gamma$ and $\delta$.

   Par.                                                                                                                             
  ------ ------ ------------------ -------------- -------------- --------------- ----------- ------------------ ------------------- ---------------
          $n$       $\hat\mu$       $\hat\sigma$   $\hat\gamma$   $\hat\delta$    $\hat\mu$     $\hat\sigma$       $\hat\gamma$      $\hat\delta$
           30         0.0052          -0.0347         0.0018         -2858.8       -0.0205        -0.0350             -0.0197           -2502.9
           50        -0.0026          -0.0269        -0.0017         -623.9        -0.0168        -0.0247             -0.0114           -476.8
          100         0.0005          -0.0096         0.0005          -5.8         -0.0047        -0.0090             -0.0037            27.0
          250         0.0010          -0.0028        -0.0001         -0.0883       0.0008         -0.0035             -0.0007           -0.0890
          500        -0.0001          -0.0003         0.0001         -0.0387       -0.0008        -0.0004             -0.0007           -0.0391
          1000   4$\times10^{-5}$     -0.0001         0.0001         -0.0202       0.0002     5$\times10^{-5}$   -5$\times10^{-5}$      -0.0187
           30         0.3303           0.0867         0.1034      7$\times10^8$    0.3054          0.0874             0.0952         5$\times10^8$
           50         0.1465           0.0493         0.0468      1$\times10^8$    0.1377          0.0490             0.0441         4$\times10^7$
          100         0.0590           0.0215         0.0191      7$\times10^4$    0.0557          0.0214             0.0181         5$\times10^6$
          250         0.0203           0.0081         0.0069         0.1494        0.0195          0.0082             0.0065            0.1505
          500         0.0102           0.0039         0.0034         0.0605        0.0096          0.0039             0.0032            0.0605
          1000        0.0050           0.0019         0.0016         0.0283        0.0047          0.0019             0.0016            0.0277
           30         0.5747           0.2965         0.3216      3$\times10^4$    0.5529          0.2977             0.3092         2$\times10^4$
           50         0.3827           0.2237         0.2165      1$\times10^4$    0.3714          0.2226             0.2103            6553.9
          100         0.2428           0.1470         0.1382          277.2        0.2360          0.1465             0.1345            2271.3
          250         0.1426           0.0902         0.0835         0.3965        0.1399          0.0906             0.0809            0.3979
          500         0.1011           0.0629         0.0583         0.2491        0.0981          0.0629             0.0566            0.2490
          1000        0.0711           0.0445         0.0410         0.1694        0.0690          0.0443             0.0398            0.1677

  : Simulation results: LTP $t$, $\delta=2$.[]{data-label="table:LTPt21"}

Applications {#NumericalExamples}
============

In this section, we present several medical applications with real data that illustrate the performance and usefulness of the proposed distributions. Throughout, we adopt the epsilon-skew parameterisation $\{a(\gamma),b(\gamma)\}=\{1-\gamma,1+\gamma\}$, $\gamma \in (-1,1)$, for the LTP distributions.

Example 1: Nerve data
---------------------

In our first example we analyse the data set reported in [@CL66] which contains $n=799$ observations rounded to the nearest half in units of 1/50 second, which correspond to the time between $800$ successive pulses along a nerve fibre. We consider two baseline distributions $s$ in (\[LTPPDF\]): a Student $t$ density with $\delta>0$ degrees of freedom (LTP $t$), and a symmetric sinh-arcsinh density [@JP09; @R15] (LTP SAS). The choice for these two baseline densities is motivated as follows. The Student-$t$ distribution is a parametric family of distributions with heavier tails than the normal ones; having the normal distribution as a limit case when $\delta \rightarrow\infty$. The behaviour of the tails of the Student-$t$ density is polynomial. On the other hand, the symmetric sinh-arcsinh density (reported in Appendix A) is a parametric density function which contains a parameter that controls the tail behaviour. This distribution can capture tails heavier or lighter than those of the normal density ($\delta\lessgtr 1$), being the normal distribution a particular case ($\delta=1$). The tails of the symmetric sinh-arcsinh density are lighter than any polynomial [@JP09]. Therefore, with these two choices of the baseline density we can cover a wide range of tail behaviours. Moreover, with the additional shape parameter $\gamma$ we also cover a wide range of shapes around the shoulders of the density. Table \[table:NerveMLE\] shows the MLEs and the Akaike information criterion (AIC) associated to these models as well as some natural competitors. We also report the estimators of the LTP Normal and the lognormal distributions, which are particular cases of the LTP SAS model. The AIC favours the LTP SAS model overall, closely followed by the LTP Normal. The MLE of the parameter $\delta$ in the LTP SAS model is larger than one, indicating that the data favour a model with lighter tails than those of the lognormal distribution. The 95% confidence intervals for the parameters $\gamma$ and $\delta$ (obtained as the 0.147-level profile likelihood intervals, see [@K85]) in the LTP SAS model are $(0.31,0.53)$ and $(1.02, 1.56)$, respectively. It is worth noticing that the confidence interval for $\delta$ only include values greater than one, which are associated to tails lighter than normal. Figures \[fig:NFM\]a–\[fig:NFM\]c show the probability plots and hazard functions corresponding to the LTP SAS, lognormal, and Weibull models, which visually illustrates the fit of these models. From Figure \[fig:NFM\]d We can observe that the fitted Gamma model produces an increasing hazard function, while the LTP SAS model produces a non monotonic hazard function with decreasing tail. This behaviour coincides with that of the fitted kernel estimation of the hazard function (which was obtained using lognormal kernels).

     Model      $\hat\mu$   $\hat\sigma$    $\hat\gamma$   $\hat\delta$         AIC
  ------------ ----------- --------------- -------------- -------------- -----------------
    LTP $t$       2.59          1.04            0.40          111.90          5401.80
    LTP SAS       2.63          1.39            0.42           1.22       [**5395.71**]{}
   LTP Normal     2.59          1.05            0.40            –             5398.45
   Log-normal     1.91          1.08             –              –             5443.70
    Weibull         –       (scale) 11.27   (shape) 1.08        –             5415.40
     Gamma          –       (scale) 9.31    (shape) 1.17        –             5411.11

  : Nerve data: Maximum likelihood estimates, AIC (best value in bold).[]{data-label="table:NerveMLE"}

  ----- -----
        
   (a)   (b)
        
   (c)   (d)
  ----- -----

Example 2: PBC data
-------------------

In this section, we analyse the popular Mayo primary biliary cirrhosis (PBC) data, reported in Appendix D from [@FH91], in order to illustrate the performance of the proposed distributions in the context of AFT models. This data set contains information about the survival time and prognostic factors for 418 patients in a study conducted at Mayo Clinic between 1974 and 1984. The survival times are reported in days together with an indicator variable associated to the status of the patient at the end of the study (0/1/2 for censored, transplant, dead). [@J03] fitted, using a semiparametric method, an AFT model with five covariates: age (in years), logarithm of the serum albumin (in mg/dl), logarithm of the serum bilirubin (in mg/dl), edema, and logarithm of the prothrombin time (in seconds). Similarly, [@D10] reports the semiparametric estimators of the AFT model with an intercept parameter as follows: $(8.692,-0.025,1.498,-0.554,-0.904,-2.822)$. We consider a maximum likelihood estimation approach of the AFT model (\[AFTModel\]) containing an intercept and LTP $t$ errors with parameters $(0, \sigma_{\varepsilon}, \gamma_{\varepsilon}, \delta_{\varepsilon})$. The estimators and the AIC values are reported in Table \[table:PBCMLE\]. We can see that the estimators obtained for the model with LTP $t$ and Log-$t$ models are close to those reported by [@D10] using a semiparametric method. The AIC values favour the models with LTP-$t$ and Log Student-$t$ (Log-$t$) errors. However, these values do not provide strong evidence to distinguish between the two models, and therefore the model choice deserves further investigation. The MLE of the skewness parameter $\gamma_{\varepsilon}$ is relatively far from zero in the Log-$t$ model. However, the inclusion of this parameter produces little effect in the estimation of the degrees of freedom $\delta_{\varepsilon}$ and the regression parameters. The 95% confidence interval for $\gamma_{\varepsilon}$ in the LTP-$t$ model is $(-0.374,0.167)$, which does not rule out the value $\gamma_{\varepsilon}=0$ as a likely value of the parameter. Then, a parsimony argument favours the model with log-Student $t$ errors (Log-$t$) in this case. Moreover, we can observe that the MLEs of $\gamma_{\varepsilon}$ in the LTP $t$ and LTP Normal model have different signs. The 95% confidence interval for $\gamma_{\varepsilon}$ in the LTP Normal model is $(-0.072,0.493)$ (which indicates that $\gamma_{\varepsilon}=0$ is an unlikely value of the parameter). The reason for this difference is that the data seem to favour a model with heavier tails than normal. The lack of flexibility in the tails and the presence of extreme observations affect the estimation of the shape parameter $\gamma_{\varepsilon}$ in the LTP Normal model by pulling out this estimator in the opposite direction. This emphasises the importance of assessing the type of flexibility required for properly modelling the data.

           Model                LTP $t$       LTP Normal       Log-$t$       Log-normal
  ------------------------ ----------------- ------------ ----------------- ------------
         Intercept               7.704          7.518           7.539          7.731
            Age                 -0.026          -0.026         -0.027          -0.025
        log(Albumin)             1.552          1.529           1.554          1.472
       log(Bilirubin)           -0.587          -0.620         -0.595          -0.606
           Edema                -0.762          -0.710         -0.706          -0.840
        log(Protime)            -2.464          -2.189         -2.313          -2.371
   $\sigma_{\varepsilon}$        0.773          0.908           0.770          0.973
   $\gamma_{\varepsilon}$       -0.133          0.190             –              –
   $\delta_{\varepsilon}$        4.446            –             5.602            –
            AIC             [**635.019**]{}    642.318     [**633.702**]{}    642.222

  : PBC data: Maximum likelihood estimates, AIC (best values in bold).[]{data-label="table:PBCMLE"}

Example 3: NCCTG Lung Cancer Data {#Lung}
---------------------------------

In this section, we revisit the popular NCCTG Lung Cancer Data. This data set contains the survival times of $n=227$ patients (the total number of patients is 228 but we have removed one patient with a missing covariate, for the sake of simplicity) with advanced lung cancer from the North Central Cancer Treatment Group. The goal of this study was to compare the descriptive information from a questionnaire applied to a group of patients against the information obtained by the patient’s physician, in terms of prognostic power [@L94]. We fit an AFT model with three covariates “age” (in years),“sex” (Male=1 Female=2), “ph.ecog” (ECOG performance score, 0=good–5=dead) as well as an intercept, with different choices for the distribution of the errors in (\[AFTModel\]). Table \[table:LungMLE\] shows the MLEs associated to each of these models together with the AIC values. The AIC favours the model with LTP logistic errors, closely followed by the model with LTP SAS errors. One explanation for this is that the estimators of the LTP SAS model indicate tails heavier than normal ($\hat\delta=0.6674$), which is a tail behaviour naturally captured by the LTP logistic distribution without additional shape parameters that control the tail. The 95% confidence intervals for $\gamma$ and $\delta$ in the LTP SAS model are $(-0.05,0.60)$ and $(0.46,0.86)$, respectively, while the corresponding confidence interval for the parameter $\gamma$ in the LTP logistic model is $(0.16, 0.62)$.

     Model     LTP SAS    LTP Normal   Log-normal     LTP logistic     Log-logistic
  ----------- ---------- ------------ ------------ ------------------ --------------
   Intercept    6.2077      6.9505       6.4940          6.5538           5.9500
      Age      -0.0068     -0.0149      -0.0191         -0.0100          -0.0082
      Sex       0.4614      0.4259       0.5219          0.4243           0.4857
    ph.ecog    -0.3824     -0.3121      -0.3557         -0.3541          -0.4042
   $\sigma$     0.4639      0.8835       1.0286          0.4847           0.5360
   $\gamma$     0.3095      0.5051         –             0.4083             –
   $\delta$     0.6674        –            –               –                –
      AIC      538.2100    545.9405     563.8323    [**536.0556**]{}     545.0486

  : NCCTG Lung Cancer data: Maximum likelihood estimates, AIC (best value in bold).[]{data-label="table:LungMLE"}

It is sometimes of interest to obtain information about the remaining life of individual cancer patients. This information is used for future planning of health care, which is of financial and medical importance. Specifically, the probability that patient $i$ survives until time $t$, given that he/she was alive at time $t_i$ is given by, $$\begin{aligned}
\label{RemainingLife}
G(t\vert t_i;\bm{\theta}) = {\mathbb P}(T\leq t\vert T> t_i) = \dfrac{G(t;\bm{\theta})-G(t_i;\bm{\theta})}{1-G(t_i;\bm{\theta})},\,\,\, t\geq t_i,\end{aligned}$$ where $G$ is the distribution under the model of interest. For an AFT model, the parameter $\bm{\theta}$ contains both the regression parameters as well as the parameters of the distribution of the errors. The simplest way to obtain an estimator of this probability consists of plugging in the MLE of $\bm{\theta}$ in (\[RemainingLife\]). The 100(1-$\alpha$)% prediction interval [@H09] for a patient that survived until time $t_i$ is $[T_i^L,T_i^R]$, which satisfies $G(T_i^L\vert t_i;\bm{\theta})=\alpha_1$ and $G(T_i^R\vert t_i;\bm{\theta})=1-\alpha_2$, with $\alpha_1 + \alpha_2=\alpha$. In our application we choose $\alpha_1=\alpha_2=0.05$, and we centre the prediction intervals at the mean of the regression model with LTP logistic errors. Figure \[fig:PredInt\] shows the 90% prediction interval for the remaining life for 10 censored patients.

We can observe from Table \[table:LungMLE\] that the estimators of the regression parameters are very similar for the different choices of the distribution of the errors. At first glance, one might think that the choice of the distribution of the residual errors has little impact on the inference. However, if the interest in on predicting the remaining life of censored patients, we may obtain different intervals for different models. For instance, Figure \[fig:Patient3\] shows how different the survival functions of the remaining life for a particular censored patient, associated to the models with LTP logistic and logistic errors, can be. This emphasises the importance of the correct specification of the distribution of the residual errors.

[c]{}

[c]{}

Discussion
==========

We have proposed a flexible class of parametric distributions (LTP) with positive support that can be used for the modelling of survival data. We have shown that some members of this class of distributions represent a flexible extension of the classical choices such as the lognormal, log-logistic, and log Student-$t$ distributions. The genesis of LTP distributions allows the user to play with different baseline log-symmetric distributions in order to properly model the tail behaviour of the data. These distributions can be used to produce models that are robust to departures from the assumption of log-symmetry. Moreover, LTP distributions preserve the ease of use of the baseline log-symmetric distribution. For instance, in models that assume lognormality, a LTP-normal can be implemented with virtually the same parsimony level. In practice, we recommend to conduct a model selection between 4-parameter LTP models and the corresponding 3- and 2-parameter submodels. Given that the parameters of LTP distributions are easily-interpreted, this model selection provides information about the features favoured by the data, such as asymmetry and tail behaviour, providing in turn more insights on the phenomenon of interest. Model selection between these nested models can be conducted either using AIC or the likelihood ratio test. The good behaviour of the MLE in this family can be established by appealing to the literature on the study of inferential properties of the family of two-piece distributions, which are linked to the proposed models via a logarithmic transformation. Confidence intervals for the model parameters can be obtained by using the profile likelihood. This approach avoids relying on asymptotic results, such as normal confidence intervals (standard errors), that may not be accurate for small or moderate sample sizes.

We conclude by pointing out possible extensions of our work. Multivariate extensions of the family of LTP distributions can be produced by using copulas. This approach has the advantage of separating the role of the parameters that control the shape of the distribution of the marginals and the dependencies between the marginals. As discussed in [@H09], the plug-in estimators considered in Section \[Lung\] may produce prediction intervals of the remaining life with a smaller coverage probability. The calibration of these intervals to improve their coverage in the context of LTP models represents an interesting research line.

Acknowledgements {#acknowledgements .unnumbered}
================

We thank the Editor and two referees for their very helpful comments. FJR gratefully acknowledges research support from EPSRC grant EP/K007521/1.

Appendix A: Some density functions {#appendix-a-some-density-functions .unnumbered}
==================================

Throughout we use the notation $t=\dfrac{x-\mu}{\sigma}$.

-   The symmetric sinh-arcsinh distribution [@JP09]: $$\begin{aligned}
    f(x;\mu,\sigma,\delta) = \dfrac{\delta}{\sigma}\phi\left[\sinh\left(\delta \operatorname{arcsinh}\left(t\right)\right)\right] \dfrac{\cosh\left(\delta \operatorname{arcsinh}\left(t\right)\right)}{\sqrt{1+t^2}},\end{aligned}$$ where $\delta>0$ controls the tails of the density, and $\phi$ is the standard normal density function. Note that for $\delta=1$, this density corresponds to the normal density.

-   The exponential power distribution: $$\begin{aligned}
    f(x;\mu,\sigma,\delta) = \dfrac{\delta}{2\sigma\Gamma(1/\delta)}\exp\left(-\vert t \vert^\delta\right).\end{aligned}$$

    where $\Gamma(\cdot)$ is the gamma function. This family contains the Laplace distribution for $\delta=1$ and the normal distribution (with variance $\sigma^2/2$) for $\delta=2$.

-   The Student-$t$ distribution: $$\begin{aligned}
    f(x;\mu,\sigma,\delta) = \dfrac{\Gamma\left(\dfrac{\delta+1}{2}\right)}{\sigma\sqrt{\pi\delta } \Gamma\left(\dfrac{\delta }{2}\right)} \left(1+\dfrac{t^2}{\delta}\right)^{-\dfrac{\delta
       +1}{2}},\end{aligned}$$ where $\Gamma(\cdot)$ is the gamma function.

-   The logistic distribution: $$\begin{aligned}
    f(x;\mu,\sigma) = \dfrac{1}{\sigma} \dfrac{\exp(-t)}{[1+\exp(-t)]^2}.\end{aligned}$$

[9]{} Azzalini, A. and Genton, M. G. (2008). Robust likelihood methods based on the skew-$t$ and related distributions. [*International Statistical Review*]{} 76: 106–129. Arellano-Valle, R. B., G[ó]{}mez, H. W. and Quintana, F. A. (2005). Statistical inference for a general class of asymmetric distributions. [*Journal of Statistical Planning and Inference*]{} 128: 427–443. Barros, M., Paula, G. A. and Leiva, V. (2008). A new class of survival regression models with heavy-tailed errors: robustness and diagnostics. [*Lifetime Data Analysis*]{} 14:316-–332. Cox D. R. and Lewis P. A. W. (1966). [*The statistical analysis of series of events*]{}. Methuem, London. Ding, Y. (2010). Some new insights about the accelerated failure time model. PhD thesis, The University of Michigan. Fern[á]{}ndez, C. and Steel, M. F. J. (1998). On Bayesian modeling of fat tails and skewness. [*Journal of the American Statistical Association*]{} 93: 359–371. Ferreira, J. T. A. S. and Steel, M. F. J. (2006). A constructive representation of univariate skewed distributions. [*Journal of the American Statistical Association*]{} 101: 823–829. Fleming, T. R. and Harrington, D. P. (1991). [*Counting Processes and Survival Analysis*]{}. Wiley: New York. Hong, Y., Meeker, W. Q. and McCalley, J. D. (2009). Prediction of the remaining life of power transformers based on left truncated and right censored lifetime data. [*Annals of Applied Statistics*]{} 2: 857-879. Jin, Z., Lin, D. Y., Wei, L. J., and Ying, Z. (2003). Rank–based inference for the accelerated failure time model. [*Biometrika*]{} 90: 341–353. Jones, M. C., and Anaya-Izquierdo K. (2010). On Parameter Orthogonality in Symmetric and Skew Models. [*Journal of Statistical Planning and Inference*]{} 141: 758–770. Jones, M. C., and Pewsey A. (2009). Sinh-arcsinh Distributions. [*Biometrika*]{} 96: 761–780. Kalbfleisch, J. G. (1985). [*Probability and Statistical Inference: Volume 2: Statistical Inference*]{} (Second Edition). Springer-Verlag, New York. Kom[á]{}rek, A. and Lesaffre, E. (2008). Bayesian accelerated failure time model with multivariate doubly interval–censored data and flexible distributional assumptions. [*Journal of the American Statistical Association*]{} 103: 523–533. Kotz, S., Kozubowski, T. J. and Podg[ó]{}rski, K. (2001). [*The Laplace Distribution and Generalizations: A Revisit with Applications to Communications, Economics, Engineering, and Finance*]{}. Birkhauser, Boston. Lawless, J. F. (2003). [*Statistical Models and Methods for Lifetime Data*]{} (2nd. Edition). John Wiley & Sons. Inc., New Jersey. Loprinzi C. L., Laurie, J. A., Wieand, H. S., Krook, J. E., Novotny, P. J., Kugler, J. W., Bartel, J., Law, M., Bateman, M., Klatt, N. E. et al. (1994). Prospective evaluation of prognostic variables from patient-completed questionnaires. [*Journal of Clinical Oncology*]{} 12: 601–607. Marchenko, Y. V. and Genton, M. G. (2010). Multivariate log-skew-elliptical distributions with applications to precipitation data. [*Environmetrics*]{} 21: 318-40. Marshall A. W. and Olkin I. (1997). A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families. [*Biometrika*]{} 84: 641–652. Meeker, W. Q. and Escobar, L. A. (1998). [*Statistical Methods for Reliability Data*]{}. Wiley, New York. Mitra, D.(2012). Likelihood inference for left truncated and right truncated censored lifetime data. PhD Thesis, McMaster University. Mitzenmacher, M. (2001). A brief history of generative models for power law and lognormal distributions. [*Internet Mathematics*]{} 1: 226–251. Mudholkar, G. S. and Hutson, A. D. (2000). The epsilon-skew-normal distribution for analyzing near-normal data. [*Journal of Statistical Planning and Inference*]{} 83: 291–309. Nadarajah, S. and Bakar, S. A. A. (2013). CompLognormal: An R Package for Composite Lognormal Distributions. [*A peer-reviewed, open-access publication of the R Foundation for Statistical Computing*]{} 97. R Core Team (2013). [*R: A language and environment for statistical computing*]{}. R Foundation for Statistical Computing, Vienna, Austria. URL <http://www.R-project.org/>. Roy, V., and Dey, D. K. (2015). Propriety of posterior distributions arising in categorical and survival models under generalized extreme value distribution. [*Statistica Sinica*]{} 24: 699–722. Rubio, F. J., Ogundimu, E. O., and Hutton, J. L. (2015). On modelling asymmetric data using two–-piece sinh–arcsinh distributions. [*Brazilian Journal of Probability and Statistics*]{}, in press. Rubio, F. J. and Steel, M. F. J. (2014). Inference in Two-Piece Location-Scale models with Jeffreys Priors (with discussion). [*Bayesian Analysis*]{} 9: 1–22. Vallejos, C. and Steel, M. F. J. (2015). Objective Bayesian Survival Analysis Using Shape Mixtures of Log-Normal Distributions. [*Journal of the American Statistical Association*]{} 110: 697–710.

[^1]: [University of Warwick, Department of Statistics, Coventry, CV4 7AL, UK.]{} E-mail: Francisco.Rubio@warwick.ac.uk

[^2]: [Virginia Tech, Department of Statistics, Blacksburg, VA 24061, USA.]{}
