The function `fmrs.gendata`

generates a data set from an FMRs model. It has the form

`fmrs.gendata(nObs, nComp, nCov, coeff, dispersion, mixProp, rho, umax, ...)`

where `n`

is the sample size, `nComp`

is the order of FMRs model, `nCov`

is the number of regression covariates, `coeff`

, `dispersion`

and `mixProp`

are the parameters of regression models, i.e. regression coefficients, dispersion (standard deviation) of the errors (sub-distributions) and mixing proportions, respectively, and `rho`

is the used in the variance-covariance matrix for simulating the design matrix `x`

, and `umax`

is the upper bound for Uniform distribution for generating censoring times.

Depending on the choice of `disFamily`

, the function `fmrs.gendata`

generates a simulated data from FMRs models. For instance, if we choose `disFamily = "norm"`

, the function ignores the censoring parameter `umax`

and generates a data set from an FMR model with Normal sub-distributions. However, if we choose `disFamily = "lnorm"`

or `disFamily = "weibull"`

, the function generates data under a finite mixture of AFT regression model with Log-Normal or Weibull sub-distributions.

The `fmrs.gendata`

function returns a list which includes a vector of responses `$y`

, a matrix of covariates `$x`

and a vector of censoring indicators `$delta`

as well as the name of sub-distribution of the mixture model.

The function `fmrs.mle`

in fmrs package provides maximum likelihood estimation for the parameters of an FMRs model. The function has the following form,

`fmrs.mle(y, x, delta, nComp, ...)`

where `y`

, `x`

and `delta`

are observations, covariates and censoring indicators respectively, and `nComp`

is the order of FMRs, `initCoeff`

, `initDispersion`

and `initmixProp`

are initial values for EM and NR algorithms, and the rest of arguments of the function are controlling parameres. The function returns a fitted FMRs model with estimates of regression parameters, standard deviations and mixing proportions of the mixture model. It also returns the log-likelihood and BIC under the estimated model, the maximum number of iterations used in EM algorithm and the type of the fitted model.

Note that one can do ridge regression by setting a value for tuning parameter of the ridge penalty other than zero in the argument `lambRidge`

.

To do the variable selection we provided the function `fmrs.varsel`

with the form

`fmrs.varsel(y, x, delta, nComp, ...)`

where `penFamily`

is the penalty including `"adplasso"`

, `"lasso"`

, `"mcp"`

, `"scad"`

, `"sica"`

and `"hard"`

, and `lambPen`

is the set of tuning parameters for components of penalty. We can run the function `fmrslme`

first and use the parameter estimates as initial values for the function `fmrs.varsel`

.

There are two approaches for specifying tuning parameters: Common and Component-wise tuning parameters. If we consider choosing common tuning parameter, we can use the BIC criteria to search through the a set of candidate values in the interval \((0,\lambda_M)\), where \(\lambda_M\) is a prespecified number. The BIC is provided by the function `fmrs.varsel`

for each candidate point and we choose the optimal \(\hat\lambda\), the one that maximizes BIC. This approach will be good for the situations with enough samples sizes. It also reduces the computational burden.

On the other hand, if we consider choosing component-wise tuning parameters we use the following function to search for optimal \((\lambda_1, \ldots, \lambda_K)\) from the set of candidate values in \((0, \lambda_M)\). The function is

`fmrs.tunsel(y, x, delta, nComp, ...)`

It is a good practice run the function `fmrs.mle`

first and use the parameter estimates as initial values in the function `fmrs.tunsel`

. The function `fmrs.mle`

then returns a set optimal tuning parameters of size `nComp`

to be used in variable selection function. Note that this approach still is under theoretical study and is not proved to give optimal values in general.

We use a simulated data set to illustrate using `fmrs`

package. We generate the covariates (design matrix) from a multivariate normal distribution of dimension `nCov=10`

and sample size 200 with mean vector \(\bf 0\) and variance-covariance matrix \(\Sigma=(0.5^{|i-j|})\). We then generate time-to-event data from a finite mixture of two components AFT regression models with Log-Normal sub-distributions. The mixing proportions are set to \(\pi=(0.3, 0.7)\). We choose \(\boldsymbol\beta_0=(2,-1)\) as the intercepts, \(\boldsymbol\beta_1=(-1, -2, 1, 2, 0 , 0, 0, 0, 0, 0)\) and \(\boldsymbol\beta_2=(1, 2, 0, 0, 0 , 0, -1, 2, -2, 3)\) as the regression coefficients in first and second component, respectively.

We start by loading necessary libraries and assigning the parameters of model as follows.

`library(fmrs)`

```
## fmrs package, Version 1.0.9, Released 2016-07-01
## Provides parameter estimation as well as variable selection in
## Finite Mixture of Accelerated Failure Time Regression and Finite
## Mixture of Regression Models.
## Furthermore, this package provides Ridge Regression and Elastic Net.
## BugReports: https://github.com/shokoohi/fmrs/issues
```

```
set.seed(1980)
nComp = 2
nCov = 10
nObs = 500
dispersion = c(1, 1)
mixProp = c(0.4, 0.6)
rho = 0.5
coeff1 = c( 2, 2, -1, -2, 1, 2, 0, 0, 0, 0, 0)
coeff2 = c(-1, -1, 1, 2, 0, 0, 0, 0, -1, 2, -2)
umax = 40
```

Using the function `fmrs.gendata`

, we generate a data set from a finite mixture of accelerated failure time regression models with above settings as follow.

```
dat <- fmrs.gendata(nObs = nObs, nComp = nComp, nCov = nCov,
coeff = c(coeff1, coeff2), dispersion = dispersion,
mixProp = mixProp, rho = rho, umax = umax,
disFamily = "lnorm")
```

Now we assume that the generated data are actually real data. We find MLE of the parameters of the assumed model using following code. Note that almost all mixture of regression models depends on initial values. Here we generate the initial values form a Normal distribution with

```
res.mle <- fmrs.mle(y = dat$y, x = dat$x, delta = dat$delta,
nComp = nComp, disFamily = "lnorm",
initCoeff = rnorm(nComp*nCov+nComp),
initDispersion = rep(1, nComp),
initmixProp = rep(1/nComp, nComp))
coefficients(res.mle)
```

```
## Comp.1 Comp.2
## Intercept -0.26344394 0.9976953
## X.1 0.32900425 -1.2803480
## X.2 0.06690871 0.1928663
## X.3 0.69669344 0.7978279
## X.4 0.51597721 0.4171002
## X.5 -0.09399042 0.1369026
## X.6 0.26679989 -0.1871274
## X.7 -0.17686523 -0.7404996
## X.8 -0.63562498 1.2610768
## X.9 1.23776252 1.1937360
## X.10 -2.20704095 -0.1349468
```

`dispersion(res.mle)`

```
## Comp.1 Comp.2
## [1,] 2.572907 3.213501
```

`mixProp(res.mle)`

```
## Comp.1 Comp.2
## [1,] 0.5693449 0.4306551
```

As we see the ML estimates of regression coefficients are not zero. Therefore MLE cannot provide a sparse solution. In order to provide a sparse solution, we use the variable selection methods developed by Shokoohi et. ale. (2016). First we need to find a good set of tuning parameters. It can be done by using component-wise tuning parameter selection function implemented in `fmrs`

as follows. In some settings, however, it is better to investigate if common tuning parameter performs better.

```
res.lam <- fmrs.tunsel(y = dat$y, x = dat$x, delta = dat$delta,
nComp = nComp, disFamily = "lnorm",
initCoeff = c(coefficients(res.mle)),
initDispersion = dispersion(res.mle),
initmixProp = mixProp(res.mle),
penFamily = "adplasso")
show(res.lam)
```

```
## An object of class 'fmrstunpar'
## Finite Mixture of Accelerated Failure Time Regression Models
## Log-Normal Sub-Distributions
## 2 Components; adplasso Penalty;
```

We have used MLE estimates as initial values for estimating the tuning parameters. Now we used the same set of values to do variable selection with adaptive lasso penalty as follows.

```
res.var <- fmrs.varsel(y = dat$y, x = dat$x, delta = dat$delta,
nComp = ncomp(res.mle), disFamily = "lnorm",
initCoeff=c(coefficients(res.mle)),
initDispersion = dispersion(res.mle),
initmixProp = mixProp(res.mle),
penFamily = "adplasso",
lambPen = slot(res.lam, "lambPen"))
coefficients(res.var)
```

```
## Comp.1 Comp.2
## Intercept -0.42767171 1.1516238
## X.1 0.28900672 -0.6005309
## X.2 0.12153691 0.1682068
## X.3 1.10559046 0.7184271
## X.4 0.44261109 0.4663552
## X.5 0.12228773 0.2349815
## X.6 0.12014270 -0.1738747
## X.7 -0.09136715 -0.4335955
## X.8 -0.68067622 0.9466941
## X.9 1.59500214 1.0218493
## X.10 -2.23858317 -0.0598671
```

`dispersion(res.var)`

```
## Comp.1 Comp.2
## [1,] 2.379193 3.181438
```

`mixProp(res.var)`

```
## Comp.1 Comp.2
## [1,] 0.5932972 0.4067028
```

The final variables that are selected using this procedure are those with non-zero coefficients. Note that a switching between components of mixture has happened here.

`round(coefficients(res.var)[-1,],3)`

```
## Comp.1 Comp.2
## X.1 0.289 -0.601
## X.2 0.122 0.168
## X.3 1.106 0.718
## X.4 0.443 0.466
## X.5 0.122 0.235
## X.6 0.120 -0.174
## X.7 -0.091 -0.434
## X.8 -0.681 0.947
## X.9 1.595 1.022
## X.10 -2.239 -0.060
```

Therefore the variable selection and parameter estimation is done simultaneously using the fmrs package.