Let X be the survival time with support X ⊂ R + and let Z = ( Z 1 , … , Z d ) ⊤ ∈ R d be a vector of covariates that may have an effect on X. Define T = min { Y , C }, where Y = log X and C is the logarithm of the censoring time, which is assumed to be noninformative and independent of X. Additionally, we use δ to denote the censoring indicator, i.e., δ = 0 if the time is censored. It is important to note that we do not directly observe X (or equivalently Y); instead, the observed data are ( T , δ ).

We say V = ( T , δ , Z ) follows a finite mixture of aft regression models of order K if the conditional density of ( T , δ ) given Z = z has the form: (2.1) f ∗ ( t , δ ; z , Ψ ) = ∑ k = 1 K π k [ f ( t ; θ k ( z ) , σ k ) ] δ [ S ( t ; θ k ( z ) , σ k ) ] 1 − δ × [ f C ( t ) ] 1 − δ [ S C ( t ) ] δ . Here, the π k s ( 0 < π k < 1, with ∑ k = 1 K π k = 1) are the mixing probabilities, and f C ( . ) and S C ( . ) are the density and survival functions of C, respectively. Note that f ( . ) and S ( . ) are the density and survival functions of Y, where θ k ( z ) = h ( β 0 k + z ⊤ β k ). In this equation, h ( . ) is a known link function, β 0 k and β k = ( β k 1 , β k 2 , … , β k d ) ⊤ are the intercept and regression coefficients, respectively, and σ k is a dispersion parameter.

It is worth noting that for each component of the mixture specified in (2.1), say the kth component, we have: Y = log X = h ( β 0 k + z ⊤ β k ) = β 0 k + z ⊤ β k + σ k ϵ . Here, ϵ has a suitable distribution such as standard normal, extreme value, generalized extreme value, or logistic. A common aft model in survival analysis is based on the Log-Normal distribution [19] in which ϵ ∼ N ( 0 , 1 ).

The vector of all parameters is: Ψ = ( β 01 , … , β 0 K , β 1 , … , β K , σ 1 , … , σ K , π 1 , … , π K − 1 ) ⊤ , which has a length of d ∗ = K ( d + 3 ) − 1, increasing with the order of the mixture.

Under the assumption of noninformative censoring, (2.2) f ∗ ( t , δ ; z , Ψ ) ∝ ∑ k = 1 K π k [ f ( t ; θ k ( z ) , σ k ) ] δ [ S ( t ; θ k ( z ) , σ k ) ] 1 − δ .

Let Y be the response variable, and let Z = ( Z 1 , … , Z d ) ⊤ be a vector of covariates that may have an effect on Y. We say that ( Y , Z ) follows an fmr model of order K [24] if the conditional density of Y given z has the form: (2.3) f ( y ; z , Ψ ) = ∑ k = 1 K π k f ( y ; θ k ( z ) , σ k ) . It should be noted that for a given sample when all the observations are failure times, i.e., there is no censoring ( δ i = 1), the density f ∗ ( t , δ ; z , Ψ ) in (2.2) is proportional to the density of a finite mixture of regression models. Therefore, (2.4) f ∗ ( y ; z , Ψ ) ∝ f ∗ ( t , δ = 1 ; z , Ψ ) . As a result, the mle of parameters of the fmr model is the same if we use the finite mixture of aft regression model with no censoring. Therefore, the fmrs package can also be used to analyze data using the fmr model.

Let τ k ( m ) = diag τ i k ( m ) : i = 1 , … , n , k = 1 , … , K. Denote the pseudo-survival times as: (3.4) t i k ( m ) = δ i t i + ( 1 − δ i ) z i ⊤ β k ( m ) + σ k ( m ) A ( ω i k ( m ) ) , where ω i k ( m ) = ( t i − z i ⊤ β k ( m ) ) / σ k ( m ) , A ( ω ) = ϕ ( ω ) / ( 1 − Φ ( ω ) ), and ϕ ( . ) and Φ ( . ) are the density and cumulative distribution functions of N ( 0 , 1 ), respectively. For each k = 1 , 2 , … , K, let T k ( m ) = ( t 1 k ( m ) , t 2 k ( m ) , … , t n k ( m ) ) ⊤ , and let Z = ( z 1 , z 2 , … , z n ) ⊤ be the n × d dimensional design matrix. The updated estimates of the parameters for the mixture of Log-Normal aft models (when t is actually the logarithm of t) and the mixture of Normal models for k = 1 , 2 , … , K are given by (3.5) β k ( m + 1 ) = Z ⊤ τ k ( m ) Z − 1 Z ⊤ τ k ( m ) T k ( m ) , and (3.6) σ k ( m + 1 ) = ∑ i = 1 n τ i k ( m ) ( t i k ( m ) − z i ⊤ β k ( m ) ) 2 ∑ i = 1 n τ i k ( m ) δ i + ( 1 − δ i ) { A ( ω i k ( m ) ) [ A ( ω i k ( m ) ) − ω i k ( m ) ] } .

There is no closed-form solution for parameter estimation in the Weibull distribution. Hence, we use a numerical method such as the Newton-Raphson algorithm [29]. The iterative algorithm is basically given as β new = β old − I − 1 ( β old ) U ( β old ) , where U and I are the first and second derivative functions of the log-likelihood, respectively, evaluated at β old .

Let Y ∗ follow a Weibull distribution. Then y = log y ∗ follows the extreme value distribution with the probability density function (pdf) and cumulative distribution function (cdf) given by f ( y ; x ⊤ β , σ ) = 1 σ exp y − x ⊤ β σ exp − exp y − x ⊤ β σ and F ( y ; x ⊤ β , σ ) = 1 − exp − exp y − x ⊤ β σ , respectively. For the kth component, we have U ( β 0 k ; β k old , y ) = − ∑ i = 1 n t i k old δ i σ k old + ∑ i = 1 n t i k old σ k old e y i − x i ⊤ β k old σ k old , and U ( β j k ; β k old , y ) = − ∑ i = 1 n t i k old δ i σ k old x i j + ∑ i = 1 n t i k old σ k old x i j e y i − x i ⊤ β k old σ k old , for j = 1 , … , d. For the second derivatives, we have I ( β j k , β r k ; β k old , y ) = − ∑ i = 1 n t i k old σ k old x i j x i r exp y i − x i ⊤ β k old σ k old and I ( β j k , β 0 k ; β k old , y ) = − ∑ i = 1 n t i k old σ k old x i j exp y i − x i ⊤ β k old σ k old , for j , r = 1 , … , d, and I ( β 0 k , β 0 k ; β k old , y ) = − ∑ i = 1 n t i k old σ k old exp y i − x i ⊤ β k old σ k old .

Algorithm 1 in Appendix A provides instructions for obtaining mles and Ridge estimators.

It is widely acknowledged that finite mixture models, specifically finite mixtures of regression models, are identifiable up to a permutation [14, 9]. As a consequence, in practical applications, it is common for the estimated components to deviate from the order of the simulation setup or the true order in the population (which remains unknown). To establish the correct order, it becomes necessary to possess some knowledge about the regression coefficients’ locations and select initial values accordingly. When analyzing real data, one possible approach to rearranging the components is to consider the order of their grand means.

Two approaches are available for choosing tuning parameters: the common approach and the component-wise approach. If the common tuning parameter approach is adopted, one can choose the value that minimizes the bic from a set of candidates within the interval ( 0 , λ max ], where λ max is a pre-specified value. This approach is suitable for datasets with sufficiently large sample sizes, and it reduces the computational burden when a common tuning parameter is used.

On the other hand, if we adopt the component-wise approach, we will search for the optimal values ( λ 1 , … , λ K ) from a set of candidate values derived from the interval ( 0 , λ max ]. We choose the combination that minimizes the component-wise bic, which is defined below.

Let τ ˜ i k be the mle in (3.2). For a given k, the component-wise log-likelihood is defined as (4.3) ℓ ˜ n k ( Ψ k ) = ∑ i = 1 n τ ˜ i k log f Y ( t i , θ k ( z i ) , σ k ) δ i S Y ( t i , θ k ( z i ) , σ k ) 1 − δ i .

Let Ψ ˆ k ( λ k ) be the mple under (4.3) for a given λ k ; i.e., Ψ ˆ k ( λ k ) = arg max Ψ k ℓ ˜ n k ( Ψ k ) − n π k ∑ j = 1 d p λ k ( | β j k | ) . We define the component-wise bic as (4.4) BIC k ( λ k ) = − 2 ℓ ˜ n k ( Ψ ˆ k ( λ k ) ) + DF ( λ k ) × log n k , where n k = ∑ i = 1 n τ ˜ i k , DF ( λ k ) = ∑ j = 1 d I ( β ˆ k j ≠ 0 ) is the number of estimated non-zero coefficients, and ℓ ˜ n k is evaluated at Ψ ˆ k ( λ k ). The component-wise tuning parameter is then chosen as λ ˆ k = arg min λ k ∈ ( 0 , λ max ] BIC k ( λ k ) , for k = 1 , … , K . By choosing a grid of values in the interval ( 0 , λ max ] with a length of L, the total number of searches required to find K tuning parameters is reduced to K × L. In contrast, using the non-component-wise approach would require L K searches to find K tuning parameters. By adopting the component-wise approach, we significantly reduce the number of searches. It is important to note that although this approach has not been theoretically studied, our simulations have shown promising results. Algorithm 2 in Appendix A provides instructions for obtaining the tuning parameters.

So far, we assumed that the order of fmrs is known a priori. However, in many real applications, such information is not available. For order selection in mixture model setups, information criteria such as bic have been extensively studied [24, 14]. In fmrs, we suggest using bic for order selection. Consider the possible values K ∈ { 1 , … , K max } for the order of mixture, where K max is a pre-specified upper bound. The optimal order, denoted as K ˆ, is chosen as the one that minimizes the quantity BIC ∗ ( K ) = − 2 ℓ n ( Ψ ˆ n , K ) + ( 3 K − 1 ) + ∑ k = 1 K ∑ j = 1 d I ( β ˆ j k ≠ 0 ) log n , where ℓ n in (3.1) is evaluated at Ψ ˆ n , K , the vector of estimated parameters with order K. Note that this approach has yet to be theoretically studied. In simulation studies, however, promising results are observed.

The command ‘help(package = “fmrs”)’ returns a list of all available arguments in fmrs. The basic functions are implemented in the C programming language, which greatly improves memory management and package speed. The package utilizes S4 objects and methods. Table 1 presents a list of all S4 classes along with their descriptions. Table 2 provides a comprehensive list of S4 methods and functions, including their arguments and associated default values. Table 3 displays a list of arguments, controls, and notations, along with their default values and descriptions. Currently, the fmaftr supports Log-Normal and Weibull distributions, while the fmr includes the Normal distribution in the package.

The fmrs package has two main classes and four main functions. The classes are as follows: •

The fmrsfit class

The package introduces 17 functions, with four of them being the main functions. The arguments for these functions are described in Tables 1-3. Below, we provide a brief introduction to these functions. •

The fmrs.gendata function

fmrs.gendata(nObs, nComp, nCov, coeff,

dispersion, mixProp, rho, umax, disFamily, ...),

where nObs, nComp, nCov, coeff, dispersion, mixProp and disFamily represent the sample size, the order of fmrs, the number of regression covariates, the regression coefficients, the dispersions of errors, the mixing proportions, and the distribution of components of fmrs, respectively.

It is important to note that rho (i.e., ρ) is used in the variance-covariance matrix to simulate the design matrix X from a multivariate Gaussian distribution with mean 0 and variance-covariance matrix Σ X = ( ρ | l − m | ). The right censoring times are generated using a Uniform distribution with lower and upper bounds of 0 and umax, respectively. Depending on the choice of disFamily, the function fmrs.gendata generates a dataset from fmaftr or fmr. The default value is disFamily = "lnorm". If disFamily = "norm", the function ignores the censoring parameter umax and generates a dataset from fmr with Normal sub-distributions. On the other hand, if disFamily = "lnorm" or disFamily = "weibull", the function produces a dataset from fmaftr with Log-Normal or Weibull sub-distribution. Consequently, fmrs.gendata returns a list containing a vector of responses y, a matrix of covariates x, a vector of censoring indicators delta, and the name of the sub-distributions of the mixture model. •

The fmrs.mle function

fmrs.mle(y, delta, x, nComp, disFamily,

initCoeff, initDispersion, initmixProp,

oracleset, ... ).

Here, delta, initCoeff, initDispersion, initmixProp, and oracleset represent the censoring indicators, the initial values for regression coefficients, the dispersions, and the mixing proportions, and the set of oracle covariates that should be included in each component of the mixture model, respectively. The remaining arguments in fmrs.mle are controlling parameters.

This function returns a fitted fmrs model that includes the mle of regression parameters, standard deviations, and mixing proportions based on the EM algorithm. The output also includes the log-likelihood and bic for the fitted model, the maximum number of iterations used in the EM algorithm, and the type of the fitted fmrs (i.e., FMAFTR or FMR). To perform Ridge regression, a positive value must be chosen for lambRidge, which is the tuning parameter of the Ridge penalty.

The default values for the arguments are as follows: nComp = 2, disFamily = "lnorm", lambRidge = 0, nIterEM = 400, nIterNR = 2, conveps = 1e-08, convepsEM = 1e-08, convepsNR = 1e-08, and NRpor = 2. •

The fmrs.tunsel function

fmrs.tunsel(y, delta, x, nComp, disFamily,

initCoeff, initDispersion, initmixProp,

penFamily, lambRidge, oracleset, lambMCP,

lambSICA, LambMin, LambMax, nLamb, ...).

Here, penFamily, lambMCP, and lambSICA represent the penalty function, the hyper-parameter for mcp, and the hyper-parameter for sica penalty. Additionally, LambMin, LambMax, and nLamb specify the minimum, maximum, and the number of tuning parameters used to obtain the optimal tuning parameter based on the component-wise tuning parameter selection approach. The function returns an fmrstunpar class.

The default values for the arguments are as follows: disFamily = "lnorm", penFamily = "lasso", lambRidge = 0, nIterEM = 400, nIterNR = 2, conveps = 1e-08, convepsEM = 1e-08, convepsNR = 1e-08, NRpor = 2, gamMixPor = 1, cutpoint = 0.05, LambMin = 0.01, LambMax = 1, and nLamb = 100. •

The fmrs.varsel function

fmrs.varsel(y, delta, x, nComp, disFamily,

initCoeff, initDispersion, initmixProp,

penFamily, lambPen, lambRidge, oracleset,

lambMCP, lambSICA, ... ).

Here, lambPen represents the set of tuning parameters for the penalty function. The function returns an fmrstfit class that stores the values of mple for the regression parameters.

The default values for the arguments are as follows: disFamily = "lnorm", penFamily = "lasso", lambRidge = 0, nIterEM = 2000, nIterNR = 2, conveps = 1e-08, convepsEM = 1e-08, convepsNR = 1e-08, NRpor = 2, gamMixPor = 1, and cutpoint = 0.05. •

Additional functions

In order to illustrate the application of fmrs, we begin by generating a dataset from an fmaft model. It is worth noting that for a comprehensive simulation study, users can refer to [29].

We generate the covariates from a multivariate normal distribution with a dimension of 10 and a sample size of 500. The mean vector is set to 0, and the variance-covariance matrix is Σ = ( 0.25 | l − m | ). Subsequently, we simulate time-to-event data from a finite mixture of two components using aft regression models with Log-Normal sub-distributions. We load the necessary libraries and assign the parameters of the model. The parameter values chosen for this simulation are provided in the following code:

One can use fmrs.gendata to generate data from an fmaftr model as follows:

Regrettably, the use of R for random generation produces distinct datasets across various machines and operating systems. Therefore, we have included alternative Python code to ensure reproducibility.

With the simulated dataset in hand, we proceed to estimate the mles of the model parameters using the fmrs.mle function. It is worth noting that the initial values for the regression parameters are generated from a standard normal distribution.

It is evident that the mles of the regression coefficients are not equal to zero. As a result, the mle approach alone cannot provide a sparse solution. To achieve sparsity, we utilize the variable selection method developed by [29]. Once we obtain the mles, the next step is to determine a set of suitable tuning parameters. This can be accomplished by employing the component-wise approach implemented in the fmrs.tunsel function. However, in certain scenarios, it is worthwhile to explore whether the common tuning parameter approach yields superior results. This can be investigated through data-driven simulation studies, for example.

We have utilized the mle estimates as initial values to obtain the tuning parameters. In this phase, the same set of values is employed to conduct variable selection with an adaptive lasso penalty.

As mentioned in Section 2, the mle and mple of an fmr model can be obtained by disregarding the censoring in the fmaftr setting. We select the following parameters to generate the data from an fmr model.

By specifying "norm" for disFamily in fmrs.gendata, we generate a dataset from an fmr model using the following code:

Similar to the above, we use Python to generate the data.

The mle of the fmr model parameters are obtained as follows:

The following code is used in selecting the component-wise tuning parameters:

Having selected the tuning parameters, we perform variable selection in fmr as follows:

We conducted a simulation study to compare the performance of fmrs with existing methods. Currently, fmrlasso is the only method available that offers a variable selection in fmr, and there is no package providing variable selection in fmaftr.

In our simulation study, we consider the model y i = x i β + e i , i = 1 , … , n, where e i ∼ i i d N ( 0 , σ 2 ). For the fmr model, we set K = 2, n = 500, d = 5, π = [ 0.4 , 0.6 ], ρ = 0.25, σ = [ 1 , 1 ], β 1 = [ − 1 , 2 , 0 , 0 , − 1 , 2 ] ⊤ , β 2 = [ 2 , 1 , 1 , 0 , 0 , 0 ] ⊤ .

We generated r = 100 datasets and evaluated the performance. The simulation results are presented in Table 4. The results demonstrate that fmrs outperforms fmrlasso in terms of correctly identifying zero or non-zero coefficients. The codes are given in the supplementary materials.

It is worth noting that fmrs significantly outperforms fmrlasso in terms of computational speed. In the aforementioned simulation study, where bic was used, fmrs was found to be over 20 times faster than fmrlasso. Moreover, unlike fmrlasso, fmrs provides variable selection capabilities for various penalties such as mcp, scad, and others. Additionally, fmrs can handle variable selection in the presence of censored observations, which is a limitation of fmrlasso.

As stated by the reviewers, the non-mixture versions of the models implemented in fmrs are highly beneficial for researchers. Therefore, we have added these models to the package. The users need to specify n C o m p = 1 to fit such models. A sample code for the accelerated failure time regression model is given as follows: