Profile and Non-Profile MM Modeling of Cluster Failure Time and Analysis of ADNI Data

Abstract

Motivated by the Alzheimer’s Disease Neuroimaging Initiative (ADNI) data, the objective of integration of important biomarkers for the early detection of Mild Cognitive Impairment (MCI) to Alzheimer’s disease (AD) as a therapeutic intervention is most likely to be beneficial in the early stages of disease progression. Developing predictors for MCI to AD comes down to genotype variables such that the dimension of predictors increases as the sample becomes large. Thus, we consider the sparsity concept of coefficients in a high-dimensional regression model with clustered failure time data such as ADNI, which enables enhancing predictive performances and facilitates the model’s interpretability. In this study, we propose two MM algorithms (profile and non-profile) for the shared frailty survival model firstly and then extend the two proposed MM algorithms to regularized estimation in sparse high-dimensional regression model. The convergence properties of our proposed estimators are also established. Furthermore simulation studies and analysis of ADNI data are illustrated by our proposed methods.

1. Introduction

In biomedical research, we often encounter clustered failure time data in which individuals from the same cluster (e.g., family) share common genetic and/or environmental factors. An illustrative example comes from the Alzheimer’s Disease Neuroimaging Initiative study where each participant may develop Mild Cognitive Impairment (MCI) and/or develop Alzheimer’s disease (AD). The times to these two events for each individual are expected to be highly associated.

In order to account for the correlation of associated failure times, shared frailty or random effect models are commonly used by researchers ([1,2,3,4]). In particular, the gamma frailty model ([5,6,7,8]) is numerically convenient for such analysis since the likelihood function exhibits a closed form. For example, ref. [2] proposed a frailty model where individuals within the same cluster share a common random effect. Refs. [9,10,11] provided a good and comprehensive review of the applications of different frailty models for clustered survival data. The authors of [12] compared several commonly used methods and demonstrated the advantages of the gamma frailty model. Refs. [13,14] also explored the nonparametric frailty model and the between-within frailty model for correlated survival data.

For high-dimension regression analysis, an important and useful strategy is to exploit sparsity and assume that the true parameters lie in a low-dimensional subspace. In the past years, sparsity-restricted estimation has attracted a great deal of attention in high-dimensional regression models. That is, only a few regression coefficients are assumed to be nonzero ([15,16]). In addition to the widely used LASSO penalty, there exists other types of regularization methods such as ridge regression ([17]), bridge regression ([18]), and elastic net ([19]). In order to improve the performance of the LASSO, many modifications such as the adaptive Lasso, the smoothly clipped absolute deviation (SCAD, [20]), and minimax concave penalty (MCP, [21]) have also been proposed.

In this paper, we develop fast and efficient algorithms for regularized estimation in the general frailty model for clustered failure time data. As the model parameters consists of both the regression coefficients and the unknown nonparametric baseline hazard function, computation in the frailty model with survival data is usually intensive, especially when frailty does not result in a closed-form likelihood function. The regularized estimation poses even more challenges for the problem and further complicates computational burden. The existing approaches of general frailty model almost rely on the EM algorithm. The EM algorithm adopts Newton’s method in its maximization step, which involves matrix inversion and may not have good performance in a high-dimensional environment. The minorization-maximization (MM) algorithm possesses an ascent property, which driving the target likelihood function to increase and reliably converges to the maximum from well-chosen initial values ([22,23,24]). In particular, ref. [25] proposed three different MM algorithms for the gamma frailty model and further demonstrated their utility for the estimation in high-dimensional situations via decomposing a high-dimensional objective function into separable low-dimensional functions. In this paper, we first develop a pair of MM algorithms, including profile and non-profile, for the general frailty model. As demonstrated by [25] in numerical studies, the decomposition inosculated with regularized estimation well in sparse high-dimensional settings. For illustration, we chose to use concave penalties such as the smoothly clipped absolute deviations penalty (SCAD, [20]) and the minimax concave penalty (MCP, [21]) to explore the sparsity because they both possess the good property of unbiasedness.

The rest of the paper is organized as follows. In Section 2, we provide the model description, an overview of MM principle, and then propose a pair of MM algorithms. In Section 3, we derive regularized estimation methods via profile and non-profile MM algorithms in sparse high-dimensional regression setting. The convergence properties of the proposed algorithms are provided in Section 4. In Section 5, a series of simulation studies was conducted to assess the finite-sample performances of the proposed methods. Section 6 provides a real application to the ADNI data.

2. The Model and Estimation

2.1. The Model

Consider datasets from some population that contains $M_i⩾1$ individuals in i-th subgroup of the population, $i=1,\dots ,B$. Individuals within the i-th subgroup have dependent event times due to some unobserved covariate information summarized in a frailty, ${\omega }_i$. Let $Y_{ij}$ be the event time, let $C_{ij}$ be the censoring time, and let ${\mathbf{X}}_{ij}^{\top }=(X_{ij1},\dots ,X_{ijp})$ denote the potential covariates for the j-th individual in the i-th subgroup. The censoring time $C_{ij}$ is assumed to be independent of $Y_{ij}$, given covarites ${\mathbf{X}}_{ij}$ and frailty ${\omega }_i$. Define $t_{ij}=min(Y_{ij},C_{ij})$ and ${\delta }_{ij}=I(Y_{ij}⩽C_{ij})$, where $I(·)$ denotes the indicator function. Suppose that censorship is noninformative. Conditional on frailty ${\omega }_i$, the hazard rate function is of the following form:

$$\begin{array}{c}\hfill \lambda (t_{ij}|{\mathbf{X}}_{ij},{\omega }_i)={\omega }_i{\lambda }_0(t_{ij})exp\left\{{\mathbf{X}}_{ij}^{\top }\mathit{\beta }\right\},\end{array}$$

(1)

where ${\lambda }_0(·)$ is an arbitrary baseline hazard rate, and $\mathit{\beta }$ is a vector of unknown parameters. Assume that ${\omega }_i,i=1,\dots ,B$ are independent and identically distributed with density function $f({\omega }_i|\mathit{\theta })$ on the domain $\mathbb{W}$. Denote $\mathit{\alpha }=(\mathit{\theta },\mathit{\beta },{\Lambda }_0)$; we then propose profile MM and non-profile MM methods to estimate parameters based on the minorization-maximization (MM) principle.

2.2. An Overview of MM Principle

Before introducing the two proposed methods, let us review the MM principle. Assume that $Y_{obs}$ is the observed data, $\ell (\mathit{\alpha }\mid Y_{obs})$ is the log-likelihood function with unknown parameter $\mathit{\alpha }={({\alpha }_1,\dots ,{\alpha }_q)}^T$, and the maximum likelihood estimate of $\mathit{\alpha }$ is $\widehat{\mathit{\alpha }}$ = argmax $\ell (\mathit{\alpha }\mid Y_{obs})$. The MM principle involves two M steps: One is a minorization step, and the other is maximization step. In a maximization problem, the first step is minorization, which aims to construct a minorization/surrogate function for the objective log-likelihood function $\ell (\mathit{\alpha }\mid Y_{obs})$ to be maximized through a series of inequalities satisfying the following two conditions:

$$\left\{\begin{array}{cc}Q(\mathit{\alpha }\mid {\mathit{\alpha }}^{(k)})\hfill & ⩽\ell (\mathit{\alpha }\mid Y_{obs}),\hfill \\ Q({\mathit{\alpha }}^{(k)}\mid {\mathit{\alpha }}^{(k)})\hfill & =\ell ({\mathit{\alpha }}^{(k)}\mid Y_{obs}),\hfill \end{array}\right.$$

(2)

where ${\mathit{\alpha }}^{(k)}$ is the k-th approximation of $\widehat{\mathit{\alpha }}$. Once the minorization function $Q(\mathit{\alpha }\mid {\mathit{\alpha }}^{(k)})$ is successfully constructed for objective function $\ell (\mathit{\alpha }\mid Y_{obs})$, the following maximization step is to maximizing the surrogate function $Q(\mathit{\alpha }\mid {\mathit{\alpha }}^{(k)})$ to obtain the $(k+1)$-th approximation of $\widehat{\mathit{\alpha }}$ rather than the objective function, i.e.,

$${\mathit{\alpha }}^{(k+1)}=\mathrm{argmax} Q(\mathit{\alpha }\mid {\mathit{\alpha }}^{(k)}).$$

By MM principle, we may have $\ell ({\mathit{\alpha }}^{(k+1)}\mid Y_{obs})⩾Q({\mathit{\alpha }}^{(k+1)}\mid {\mathit{\alpha }}^{(k)})⩾Q({\mathit{\alpha }}^{(k)}\mid {\mathit{\alpha }}^{(k)})=\ell ({\mathit{\alpha }}^{(k)}\mid Y_{obs})$, and the values of the objective function continue to increase until convergence.

2.3. Profile MM Estimation Procedure

For general shared frailty models, we can write the observed log-likelihood function $\ell (\mathit{\alpha }|Y_{obs})$ as follows:

$$\begin{array}{c}\hfill \ell (\mathit{\alpha }|Y_{obs})=\sum _{i=1}^{B}log{\int }_{\mathbb{W}}{\tau }_i({\omega }_i|\mathit{\alpha })\mathrm{d}{\omega }_i,\end{array}$$

(3)

where the following is the case.

$${\tau }_i({\omega }_i|\mathit{\alpha })=f({\omega }_i|\mathit{\theta })\prod _{j=1}^{M_i}\left\{{\left[{\lambda }_0(t_{ij}){\omega }_{i}exp({\mathbf{X}}_{ij}^{\top }\mathit{\beta })\right]}^{{\delta }_{ij}}exp\left[-{\Lambda }_0(t_{ij}){\omega }_{i}exp({\mathbf{X}}_{ij}^{\top }\mathit{\beta })\right]\right\}.$$

Generally, the Laplace transform of the frailty’s distribution is theoretically intractable. Hence, the explicit form of marginal hazard is not available to us. However, that has not stopped us from developing the profile MM approach for estimation in general shared frailty model. Define the following:

$$v_i({\omega }_i|{\mathit{\alpha }}^{(k)})=\frac{{\tau }_i({\omega }_i|{\mathit{\alpha }}^{(k)})}{{\int }_{\mathbb{W}}{\tau }_i({\omega }_i|{\mathit{\alpha }}^{(k)})\mathrm{d}{\omega }_i},$$

and rewrite the objective function as follows.

$$\begin{array}{c}\hfill \ell (\mathit{\alpha }|Y_{obs})=\sum _{i=1}^{B}log\left[{\int }_{\mathbb{W}}\frac{{\tau }_i({\omega }_i|\mathit{\alpha })}{v_i({\omega }_i|{\mathit{\alpha }}^{(k)})}·v_i({\omega }_i|{\mathit{\alpha }}^{(k)})\mathrm{d}{\omega }_i\right].\end{array}$$

(4)

By Jensen’s inequality, we have the following:

$$\phi \left[{\int }_{\mathbb{X}}h(x)·g(x)\mathrm{d}x\right]⩾{\int }_{\mathbb{X}}\phi \left[h(x)\right]·g(x)\mathrm{d}x,$$

where $\mathbb{X}$ is a subset of the real line $\mathbb{R}$, $\phi (·)$ is the concave function, $h(·)$ is an arbitrary real-valued function defined on $\mathbb{X}$, and $g(·)$ is a density function defined on $\mathbb{X}$. In Equation (4), $v_i({\omega }_i|{\mathit{\alpha }}^{(k)})$ is a density function, and choosing $h(x)$ as ${\tau }_i({\omega }_i|\mathit{\alpha })/v_i({\omega }_i|{\mathit{\alpha }}^{(k)})$, we can apply the above Jensen’s inequality and construct the following surrogate function for $\ell (\mathit{\alpha }|Y_{obs})$:

$$\begin{array}{c}\hfill Q_1(\mathit{\alpha }|{\mathit{\alpha }}^{(k)})=Q_{11}(\mathit{\theta }|{\mathit{\alpha }}^{(k)})+Q_{12}(\mathit{\beta },{\Lambda }_0|{\mathit{\alpha }}^{(k)}),\end{array}$$

(5)

where the following is the case:

$$\begin{array}{c}\hfill Q_{11}(\mathit{\theta }|{\mathit{\alpha }}^{(k)})=\sum _{i=1}^B{\int }_{\mathbb{W}}log\left[f({\omega }_i|\mathit{\theta })\right]·v_i({\omega }_i|{\mathit{\theta }}^{(k)},{\mathit{\alpha }}^{(k)})\mathrm{d}{\omega }_i,\end{array}$$

(6)

which only consists of parameter $\mathit{\theta }$. Moreover, we have the following:

$$\begin{array}{c}\hfill Q_{12}(\mathit{\beta },{\Lambda }_0|{\mathit{\alpha }}^{(k)})=\sum _{i=1}^B\sum _{j=1}^{M_i}\left[{\delta }_{ij}log({\lambda }_0(t_{ij}))+{\delta }_{ij}{\mathbf{X}}_{ij}^{\top }\mathit{\beta }-A_i^{(k)}{\Lambda }_0(t_{ij})exp({\mathbf{X}}_{ij}^{\top }\mathit{\beta })\right],\end{array}$$

(7)

where $A_i^{(k)}={\int }_{\mathbb{W}}{\omega }_i·v_i({\omega }_i|{\mathit{\alpha }}^{(k)})\mathrm{d}{\omega }_i,i=1,\dots ,B.$ The minorizing function $Q_1(\mathit{\alpha }|{\mathit{\alpha }}^{(k)})$ separates parameters $\mathit{\theta }$ and $(\mathit{\beta },{\Lambda }_0)$ into (6) and (7), respectively. In the second M-step, the updated estimates of $\mathit{\theta }$ involve maximizing (6) numerically. Due to the presence of the nonparametric component ${\Lambda }_0$, updating $(\mathit{\beta },{\Lambda }_0)$ is still a big challenge. Following [7], we apply the profile estimation method to profile out ${\Lambda }_0$ in $Q_{12}(\mathit{\beta },{\Lambda }_0|{\mathit{\alpha }}^{(k)})$, which results in the estimate of ${\Lambda }_0$ given $\mathit{\beta }$.

$$\begin{array}{c}\hfill {\displaystyle \mathrm{d}{\widehat{\Lambda }}_0(t_{ij})=\frac{{\delta }_{ij}}{{\sum }_{r=1}^B{\sum }_{s=1}^{M_r}I(t_{rs}⩾t_{ij})A_r^{(k)}exp({\mathbf{X}}_{rs}^{\top }\mathit{\beta })}.}\end{array}$$

(8)

Substituting (8) into $Q_{12}(\mathit{\beta },{\Lambda }_0|{\mathit{\alpha }}^{(k)})$, we obtain the following:

$$\begin{array}{c}\hfill Q_{13}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})=\sum _{i=1}^B\sum _{j=1}^{M_i}\left\{{\delta }_{ij}{\mathbf{X}}_{ij}^{\top }\mathit{\beta }-{\delta }_{ij}log\left[\sum _{r=1}^B\sum _{s=1}^{M_r}I(t_{rs}⩾t_{ij})A_r^{(k)}exp({\mathbf{X}}_{rs}^{\top }\mathit{\beta })\right]\right\},\end{array}$$

(9)

which involves $\mathit{\beta }$ only. In Equation (9), Newton’s method and large matrix inversion are required to update $\mathit{\beta }$ in this procedure when there exist a large number of covariates. Here, we further construct minorizing functions for $Q_{13}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ to separate the regression parameters ${\beta }_1,\dots ,{\beta }_q$ from each other under the MM principle. We first use the supporting hyperplane inequality:

$$-log(x)⩾-log(x_0)-\frac{x-x_0}{x_0}$$

to minorize $Q_{13}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})$; then, we have the following surrogate function:

$$\begin{array}{c}\hfill Q_{14}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})=\sum _{i=1}^B\sum _{j=1}^{M_i}\left[{\delta }_{ij}{\mathbf{X}}_{ij}^{\top }\mathit{\beta }-\frac{{\delta }_{ij}{\sum }_{r=1}^B{\sum }_{s=1}^{M_r}I(t_{rs}⩾t_{ij})A_r^{(k)}exp({\mathbf{X}}_{rs}^{\top }\mathit{\beta })}{{\sum }_{r=1}^B{\sum }_{s=1}^{M_r}I(t_{rs}⩾t_{ij})A_r^{(k)}exp({\mathbf{X}}_{rs}^{\top }{\mathit{\beta }}^{(k)})}\right]+c,\end{array}$$

where c is a constant not depending on $\mathit{\beta }$. We apply Jensen’s inequality to the concave function $-exp(·)$ in $Q_{14}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ by rewriting the following:

$${\mathbf{X}}_{rs}^{\top }\mathit{\beta }=\sum _{p=1}^q{\pi }_{prs}[{\pi }_{prs}^{-1}{X}_{prs}({\beta }_p-{\beta }_p^{(k)})+{\mathbf{X}}_{rs}^{\top }{\mathit{\beta }}^{(k)}],$$

where ${\pi }_{prs}=|X_{prs}|/{\sum }_{p=1}^q\left|X_{prs}\right|$. Finally, the minorizing function for $Q_{14}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ is as follows:

$$\begin{array}{c}\hfill Q_{15}({\beta }_1,\dots ,{\beta }_q|{\mathit{\alpha }}^{(k)})\widehat{=}\sum _{p=1}^q{Q}_{15p}({\beta }_p|{\mathit{\alpha }}^{(k)}),\end{array}$$

(10)

where the following is the case.

$$\begin{array}{ccc}\hfill & & Q_{15p}({\beta }_p|{\mathit{\alpha }}^{(k)})=\sum _{i=1}^B\sum _{j=1}^{M_i}\{{\delta }_{ij}{X}_{pij}{\beta }_p\hfill \\ & & -\frac{{\delta }_{ij}{\sum }_{r=1}^B{\sum }_{s=1}^{M_r}I(t_{rs}⩾t_{ij})A_r^{(k)}{\pi }_{prs}exp\left[{\pi }_{prs}^{-1}{X}_{prs}({\beta }_p-{\beta }_p^{(k)})+{\mathbf{X}}_{rs}^{\top }{\mathit{\beta }}^{(k)}\right]}{{\sum }_{r=1}^B{\sum }_{s=1}^{M_r}I(t_{rs}⩾t_{ij})A_r^{(k)}exp({\mathbf{X}}_{rs}^{\top }{\mathit{\beta }}^{(k)})}\},p=1,\dots ,q.\hfill \end{array}$$

(11)

In general, the structure of minorizing function for the objective log-likelihood function $\ell (\mathit{\alpha }|Y_{obs})$ in (3) is the following:

$$\begin{array}{c}\hfill Q_{pro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})=Q_{11}(\mathit{\theta }|{\mathit{\alpha }}^{(k)})+\sum _{p=1}^q{Q}_{15p}({\beta }_p|{\mathit{\alpha }}^{(k)}),\end{array}$$

(12)

with an explicit form update of $\mathrm{d}{\Lambda }_0$ by (8). We may observe that the objective function to be maximized is decomposed into a sum of $q+1$ univariate functions from (12) as $Q_{11}(\mathit{\theta }|{\mathit{\alpha }}^{(k)})$ only consists of one parameter usually. The next maximization step of this MM algorithm only involves $q+1$ separate univariate optimizations and matrix inversion is unnecessary. Note that the success of the above profile MM algorithm requires the convergence of two improper integrals such as ${\int }_{\mathbb{W}}{\omega }_i·v_i({\omega }_i|{\mathit{\alpha }}^{(k)})\mathrm{d}{\omega }_i$ and ${\int }_{\mathbb{W}}log\left[f({\omega }_i|\mathit{\theta })\right]·v_i({\omega }_i|{\mathit{\alpha }}^{(k)})\mathrm{d}{\omega }_i$. Moreoveer, the convergences of these improper integrals are obvious when we assume the distribution of random effect comes from an exponential distribution family. The estimation proceeds by profile MM algorithm are summarized as follows:

Step 1.

Given initial values for $\mathit{\theta }$, $\mathit{\beta }$, and ${\Lambda }_0$;

Step 2.

Update $\mathit{\theta }$ via (6). Update each ${\beta }_p$ via (11) for $p=1,\dots ,q$;

Step 3.

Based on the update of $\mathit{\beta }$, compute the estimates of ${\Lambda }_0(t_{ij})$ via (8);

Step 4.

Iterate steps 2 and 3 until convergence.

2.4. Non-Profile MM Estimation Procedure

In this subsection, we bypass the profile estimation procedure in previous subsection and continue to develop new MM procedures for Equation (7) to separate parameters $\mathit{\beta }$ and nuisance baseline hazard rate ${\Lambda }_0$. Actually, to separate $\mathit{\beta }$ and ${\Lambda }_0$ of (7) is to deal with the last term $-{\Lambda }_0(t_{ij})exp({\mathbf{X}}_{ij}^{\top }\mathit{\beta })$. As in [26], we use the following arithmetic-geometric mean inequality.

$$\begin{array}{c}\hfill -\prod _{i=1}^n{x}_i^{a_i}⩾-\sum _{i=1}^n\frac{a_i}{{\left|\right|\mathbf{a}\left|\right|}_1}{x}_i^{{\left|\right|\mathbf{a}\left|\right|}_1}\end{array}$$

(13)

Here, $x_i$ and $a_i$ are non-negative. Choosing $x_1={\Lambda }_0(t_{ij})/{\Lambda }_0^{(k)}(t_{ij})$ and $x_2=exp({\mathbf{X}}_{ij}^{\top }\mathit{\beta })/$$exp({\mathbf{X}}_{ij}^{\top }{\mathit{\beta }}^{(k)})$ in inequality (13), we obtain the following surrogate function for (7):

$$\begin{array}{ccc}\hfill & & Q_2(\mathit{\beta },{\Lambda }_0|{\mathit{\alpha }}^{(k)})\hfill \\ \hfill & =& \sum _{i=1}^B\sum _{j=1}^{M_i}\left[{\delta }_{ij}log({\lambda }_0(t_{ij}))+{\delta }_{ij}{\mathbf{X}}_{ij}^{\top }\mathit{\beta }-\frac{A_i^{(k)}exp({\mathbf{X}}_{ij}^{\top }{\mathit{\beta }}^{(k)})}{2{\Lambda }_0^{(k)}(t_{ij})}{\Lambda }_0{(t_{ij})}^2-\frac{A_i^{(k)}{\Lambda }_0^{(k)}(t_{ij})}{2exp({\mathbf{X}}_{ij}^{\top }{\mathit{\beta }}^{(k)})}exp(2{\mathbf{X}}_{ij}^{\top }\mathit{\beta })\right],\hfill \\ & \widehat{=}& Q_{21}({\Lambda }_0|{\mathit{\alpha }}^{(k)})+Q_{22}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})\hfill \end{array}$$

(14)

where the following is the case.

$$\begin{array}{c}\hfill Q_{21}({\Lambda }_0|{\mathit{\alpha }}^{(k)})=\sum _{i=1}^B\sum _{j=1}^{M_i}\left[{\delta }_{ij}log({\lambda }_0(t_{ij}))-\frac{A_i^{(k)}exp({\mathbf{X}}_{ij}^{\top }{\mathit{\beta }}^{(k)})}{2{\Lambda }_0^{(k)}(t_{ij})}{\Lambda }_0{(t_{ij})}^2\right],\end{array}$$

(15)

$$\begin{array}{c}\hfill Q_{22}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})=\sum _{i=1}^B\sum _{j=1}^{M_i}\left[{\delta }_{ij}{\mathbf{X}}_{ij}^{\top }\mathit{\beta }-\frac{A_i^{(k)}{\Lambda }_0^{(k)}(t_{ij})}{2exp({\mathbf{X}}_{ij}^{\top }{\mathit{\beta }}^{(k)})}exp(2{\mathbf{X}}_{ij}^{\top }\mathit{\beta })\right].\end{array}$$

(16)

In order to obtain the nonparametric estimate of ${\Lambda }_0$, the maximization of (15) is required. For ease of computation, a one-step late skill is applied to the first order derivative of Equation (15); then, we obtain the estimate of ${\lambda }_0$ by the following:

$$\begin{array}{c}\hfill {\displaystyle \mathrm{d}{\widehat{\Lambda }}_0(t_{ij})=\frac{{\delta }_{ij}}{{\sum }_{r=1}^B{\sum }_{s=1}^{M_r}I(t_{rs}⩾t_{ij})A_r^{(k)}exp({\mathbf{X}}_{rs}^{\top }{\mathit{\beta }}^{(k)})},}\end{array}$$

(17)

which is same as (8) in the profile estimation method. To update $\mathit{\beta }$, a similar technique as dealing with $Q_{14}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ is used. We apply Jensen’s inequality to the concave function $-exp(·)$ in $Q_{22}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ by rewriting the following:

$$2{\mathbf{X}}_{ij}^{\top }\mathit{\beta }=\sum _{p=1}^q{\pi }_{pij}\left[2{\pi }_{pij}^{-1}{X}_{pij}({\beta }_p-{\beta }_p^{(k)})+2{\mathbf{X}}_{ij}^{\top }{\mathit{\beta }}^{(k)}\right],$$

where ${\pi }_{pij}=|X_{pij}|/{\sum }_{p=1}^q\left|X_{pij}\right|$. In the end, the minorizing function for $Q_{22}(\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ is as follows:

$$\begin{array}{c}\hfill Q_{23}({\beta }_1,\dots ,{\beta }_q|{\mathit{\alpha }}^{(k)})\widehat{=}\sum _{p=1}^q{Q}_{23p}({\beta }_p|{\mathit{\alpha }}^{(k)}),\end{array}$$

(18)

where the following is obtained:

$$\begin{array}{c}\hfill Q_{23p}({\beta }_p|{\mathit{\alpha }}^{(k)})=\sum _{i=1}^B\sum _{j=1}^{M_i}\left\{{\delta }_{ij}{X}_{pij}{\beta }_p-\frac{{\pi }_{pij}{A}_i^{(k)}{\Lambda }_0^{(k)}(t_{ij})exp\left[2{\pi }_{pij}^{-1}{X}_{pij}({\beta }_p-{\beta }_p^{(k)})+2{\mathbf{X}}_{ij}^{\top }{\mathit{\beta }}^{(k)}\right]}{2exp({\mathbf{X}}_{ij}^{\top }{\mathit{\beta }}^{(k)})}\right\},\end{array}$$

(19)

for $p=1,\dots ,q.$ As a result, we construct the surrogate function for the objective log-likelihood function via a non-profile MM principle as follows:

$$\begin{array}{c}\hfill Q_{nonpro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})=Q_{11}(\mathit{\theta }|{\mathit{\alpha }}^{(k)})+\sum _{p=1}^q{Q}_{23p}({\beta }_p|{\mathit{\alpha }}^{(k)}),\end{array}$$

(20)

with explicit form update of $\mathrm{d}{\Lambda }_0$ by (17). From (20), we can find similar nice features as $Q_{pro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})$; that is, $Q_{nonpro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ is a sum of $q+1$ univariate functions, which means that the next maximization (second M) step only involves $q+1$ simple univariate optimizations. It is worth noting that the parameter separated feature in the proposed profile MM and non-profile MM algorithms will help incoporate with the existing simple off-the-shelf accelerators well and brings about great effectiveness in computation time, as discussed in [25]. The estimation proceeds by non-profile MM algorithm are summarized as follows:

Step 1.

Given initial values of $\mathit{\theta }$, $\mathit{\beta }$, and ${\Lambda }_0$;

Step 2.

Update the estimate of $\mathit{\theta }$ via (6). Update the estimate of ${\beta }_p$ based on (19) for $p=1,\dots ,q$;

Step 3.

Using the updated estimate of $\mathit{\beta }$, compute the estimates of ${\Lambda }_0(t_{ij})$ via (17);

Step 4.

Iterate steps 2 and 3 until convergence.

3. Regularized Estimation Methods via MM Methods

Followed by Section 3, both proposed estimation approaches created nice parameter-separated surrogate functions to be maximized which may incorporate with regularization cleverly. In this part, we propose regularized estimation approaches based on MM algorithms in regression analysis under general shared frailty model. Many variable selection criteria arise as special cases of the general formulation as discussed in [20], where the penalized likelihood function takes the form

$$\begin{array}{c}\hfill {\ell }^P(\mathit{\alpha }|Y_{obs})=\ell (\mathit{\theta },\mathit{\beta },{\Lambda }_0|Y_{obs})-N\sum _{p=1}^{q}P(|{\beta }_p|,\lambda ),\end{array}$$

(21)

where $\ell (\mathit{\theta },\mathit{\beta },{\Lambda }_0|Y_{obs})$ is the log-likelihood function for the shared frailty model, q is the dimension of $\mathit{\beta }$, $P(·,\lambda )$ is a given nonnegative penalty function, and $\lambda ⩾0$ is a tuning parameter, which is allowed to use ${\lambda }_p$ in more general cases. This penalty shrinks some of the coefficients to zero. Under the scope of general frailty topic, the computation of MLEs is extremely complicated and hard in terms of accuracy as parameters of interests involve three blocks $\mathit{\theta }, \mathit{\beta },$ and ${\Lambda }_0$, and even more complicated when there exists a large numbers of coefficients. It is worth mentioning that the proposed profile and non-profile MM algorithms decomposed the coefficient vector $\mathit{\beta }$ from the other two blocks $\mathit{\theta },$ and ${\Lambda }_0$ and different coefficient parameters are separated from each other. This nice feature of the proposed profile and non-profile MM algorithms may mesh well with various regularization problems in (21) to produce more sparse and accurate estimates. Under profile MM algorithmic technique in (12), we obtain the corresponding minorization function for (21) as follows.

$$\begin{array}{c}\hfill Q_{pro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})-N\sum _{p=1}^{q}P(|{\beta }_p|,\lambda ).\end{array}$$

(22)

Under the non-profile MM algorithmic technique in (20), the minorization function for (21) is as follows.

$$\begin{array}{c}\hfill Q_{nonpro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})-N\sum _{p=1}^{q}P(|{\beta }_p|,\lambda ).\end{array}$$

(23)

When $P(·,\lambda )$ is piecewise differentiable, nondecreasing, and concave on $(0,\infty )$ such as $L_1$, MCP and SCAD penalties [27], explored a connection between the local quadratic approximation with MM algorithm. The penalty term $-P(·,\lambda )$ can be minorized by a local quadratic approximation form as follows:

$$\begin{array}{c}\hfill -P(|{\beta }_p|,\lambda )⩾-P(|{\beta }_p^{(k)}|,\lambda )-\frac{[{\beta }_p^2-{({\beta }_p^{(k)})}^2]P^{\prime }(|{\beta }_p^{(k)}{|}_{+},\lambda )}{2|{\beta }_p^{(k)}|}\widehat{=}\Phi ({\beta }_p|{\beta }_p^{(k)}),\end{array}$$

(24)

which is actually a one-step minorizing process. By combining function (24) with (12) or (20), respectively, we obtain the final surrogate functions for penalized likelihood function (21) as follows:

$$\begin{array}{c}\hfill Q_{pro}^P(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})=Q_{pro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})-N\sum _{p=1}^q\frac{{\beta }_p^2·P^{\prime }(|{\beta }_p^{(k)}{|}_{+},\lambda )}{2|{\beta }_p^{(k)}|}+c_1,\end{array}$$

(25)

and the following is obtained.

$$\begin{array}{c}\hfill Q_{nonpro}^P(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})=Q_{nonpro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})-N\sum _{p=1}^q\frac{{\beta }_p^2·P^{\prime }(|{\beta }_p^{(k)}{|}_{+},\lambda )}{2|{\beta }_p^{(k)}|}+c_2.\end{array}$$

(26)

Both Equations (25) and (26) are written as a sum of a series of univariate functions so that maximizing (25) and (26) will be easier than directly maximizing (21). Moreover, some simple off-the-shelf accelerators may also be used here to make the optimization problems more simplier and efficient. Regularized estimation proceeds by (profile/non-profile) MM algorithms and are summarized as follows:

Step 1.

Given initial values of $\mathit{\theta }$, $\mathit{\beta }$ and ${\Lambda }_0$;

Step 2.

Update the estimate of $\mathit{\theta }$ via (6);

Step 3.

For profile MM method, update $\mathit{\beta }$ by maximizing

$$\sum _{p=1}^q{Q}_{15p}({\beta }_p|{\mathit{\alpha }}^{(k)})-N\sum _{p=1}^q\frac{{\beta }_p^2·P^{\prime }(|{\beta }_p^{(k)}{|}_{+},\lambda )}{2|{\beta }_p^{(k)}|}.$$

For non-profile MM method, update $\mathit{\beta }$ by maximizing

$$\sum _{p=1}^q{Q}_{23p}({\beta }_p|{\mathit{\alpha }}^{(k)})-N\sum _{p=1}^q\frac{{\beta }_p^2·P^{\prime }(|{\beta }_p^{(k)}{|}_{+},\lambda )}{2|{\beta }_p^{(k)}|};$$

Step 4.

Using the updated estimate of $\mathit{\beta }$ in Step 3, compute the estimates of ${\Lambda }_0(t_{ij})$ via (8) for profile MM method and via (17) for non-profile MM method, respectively;

Step 5.

Iterate steps 2 to 4 until convergence.

Model Selection

From the recent literature, tuning parameter $\lambda$ may be selected by multiple data driven model selection criteria, such as Bayesian information criterion BIC ([28]) and generalized cross-validation GCV ([29]). In this paper, we consider a widely used BIC-type criterion, defined by the following:

$$\begin{array}{c}\hfill {BIC}_{\lambda }=-2\ell (\widehat{\mathit{\alpha }})+C_n(\widehat{S}+1)log(N),\end{array}$$

(27)

to select the tuning parameter $\lambda$, where $C_n=max\{1,log[log(q+1)]\}$, q is the dimension of $\mathit{\beta }$, and the degrees of freedom $\widehat{S}$ are defined as the number of nonzero parameters in $\widehat{\mathit{\beta }}$.

4. Theoretical Properties

We first consider the convergence properties for the profile and non-profile MM algorithms of maximizing $\ell (\mathit{\alpha })$ and ${\ell }^P(\mathit{\alpha })$. As discussed in Section 3 and Section 4, $Q_{pro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ and $Q_{nonpro}(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ are the constructed minorizing functions of $\ell (\mathit{\alpha })$ via the profile and non-profile MM approach, $Q_{pro}^P(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ and $Q_{nonpro}^P(\mathit{\theta },\mathit{\beta }|{\mathit{\alpha }}^{(k)})$ are the constructed minorizing functions of penalized likelihood function ${\ell }^P(\mathit{\alpha })$ via the profile and non-profile MM approach, $\mathit{\theta }$ and $\mathit{\beta }$ are the parameter vectors, and ${\mathit{\alpha }}^{(k)}$ is its current estimate. In this paper, we assume that the common component $Q_{11}(\mathit{\theta }|{\mathit{\alpha }}^{(k)})$ of all minorizing functions is strictly concave with respect to $\mathit{\theta }$, and the second M step of MM principle is based on the Newton–Raphson method. As strict concavity holds for all surrogate functions $Q_{pro}(·|{\mathit{\alpha }}^{(k)})$, $Q_{nonpro}(·|{\mathit{\alpha }}^{(k)})$, $Q_{pro}^P(·|{\mathit{\alpha }}^{(k)})$ and $Q_{nonpro}^P(·|{\mathit{\alpha }}^{(k)})$ by using concave penalties $L_1$, MCP and SCAD, we denote their unique maximizers by $M_{pro}({\mathit{\alpha }}^{(k)})$, $M_{nonpro}({\mathit{\alpha }}^{(k)})$, $M_{pro}^P({\mathit{\alpha }}^{(k)})$, and $M_{nonpro}^P({\mathit{\alpha }}^{(k)})$, respectively. Following Proposition 15.4.3 of Lange [30], we provide the following convergence properties.

Proposition 1.

Assume the differentiability and coerciveness of $-\ell (\mathit{\alpha })$ hold, all stationary points of $-\ell (\mathit{\alpha })$ are isolated and the subsets $\{\mathit{\alpha }\in \Omega :\ell (\mathit{\alpha })⩾\ell ({\mathit{\alpha }}^{(k)})\}$ of parameter domain Ω are compact. Then, the profile iteration sequence $({\mathit{\theta }}^{(k+1)},{\mathit{\beta }}^{(k+1)})=M_{pro}({\mathit{\alpha }}^{(k)})$ together with ${\Lambda }_0^{(k+1)}$ in (8) and non-profile sequence $({\mathit{\theta }}^{(k+1)},{\mathit{\beta }}^{(k+1)})=M_{nonpro}({\mathit{\alpha }}^{(k)})$ together with ${\Lambda }_0^{(k+1)}$ in (17) converge to the stationary point of $\ell (\mathit{\alpha })$. If the strict concavity of $\ell (\mathit{\alpha })$ also hold, then the profile MM sequence and non-profile MM sequence converge to the same maximum point of $\ell (\mathit{\alpha })$.

Proposition 2.

Suppose $P(·|\lambda )$ is piecewise differentiable, nondecreasing and concave on $(0,\infty )$, continuous at 0 and $P^{\prime }(0_{+}|\lambda )<\infty$. Then, for all ${\beta }_p\ne 0$, $\Phi ({\beta }_q|{\beta }_q^{(k)})$ as defined in (24) minorizes $-P(|{\beta }_q\left|\right|\lambda )$ at points $\pm |{\beta }_q|$. In particular, conditions in construction minorization function hold, that is, $-P(|{\beta }_q\left|\right|\lambda )⩾\Phi ({\beta }_q|{\beta }_q^{(k)})$ for all ${\beta }_q$ with equality at ${\beta }_q={\beta }_q^{(k)}$ and the ascent property $-P(|{\beta }_q^{(k+1)}\left|\right|\lambda )⩾\Phi ({\beta }_q^{(k+1)}|{\beta }_q^{(k)})⩾\Phi ({\beta }_q^{(k)}|{\beta }_q^{(k)})=-P(|{\beta }_q^{(k)}\left|\right|\lambda )$ hold.

Based on Proposition 2, we can easily obtain the following convergence properties of profile MM and non-profile MM algorithms for maximizing penalized likelihood function ${\ell }^P(\mathit{\alpha })$ defined in (21).

Proposition 3.

Assume the differentiability and coerciveness of $-{\ell }^P(\mathit{\alpha })$ hold, all stationary points of $-{\ell }^P(\mathit{\alpha })$ are isolated and the subsets $\{\mathit{\alpha }\in \Omega :{\ell }^P(\mathit{\alpha })⩾{\ell }^P({\mathit{\alpha }}^{(k)})\}$ of parameter domain Ω are compact. Then, the profile iteration sequence $({\mathit{\theta }}^{(k+1)},{\mathit{\beta }}^{(k+1)})=M_{pro}^P({\mathit{\alpha }}^{(k)})$ together with ${\Lambda }_0^{(k+1)}$ in (8) and non-profile sequence $({\mathit{\theta }}^{(k+1)},{\mathit{\beta }}^{(k+1)})=M_{nonpro}^P({\mathit{\alpha }}^{(k)})$ together with ${\Lambda }_0^{(k+1)}$ in (17) converge to the stationary point of ${\ell }^P(\mathit{\alpha })$. If the strict concavity of ${\ell }^P(\mathit{\alpha })$ also hold, then the profile MM sequence and non-profile MM sequence converge to the same maximum point of ${\ell }^P(\mathit{\alpha })$.

5. Numerical Examples

Example 1.

We independently simulate data from three different frailty models:

$$\begin{array}{c}\hfill \lambda (t|{\mathbf{X}}_{ij},{\omega }_i)={\omega }_i{\lambda }_0(t)exp\left\{{\mathbf{X}}_{ij}^{\top }\mathit{\beta }\right\},{\omega }_i\sim \left\{\begin{array}{cc}Log-normal(0,\theta ),\hfill & \theta =0.25,\hfill \\ InverseGaussian(\theta ,{\theta }^2),\hfill & \theta =1,\hfill \\ Gamma(1/\theta ,1/\theta ),\hfill & \theta =2,\hfill \end{array}\right.\end{array}$$

with three different sample sizes $(B,M)=\{(15,20),(30,13),(50,10)\}$. The true value of regression vector $\mathit{\beta }$ is set to be $(-2_6^{\top },-1_6^{\top },1_6^{\top },2_6^{\top },3_6^{\top })$ with dimension $q=30$ and all $X_i$’s are generated from independent uniform distribution between 0 and $0.5$. The censoring times are generated from independent uniform distribution to yield censoring proportion at around 15% or 30%. In this example, we numerically illustrate the efficiency of two proposed profile and non-profile MM algorithms under three different (Log-normal, Inverse Gaussian, and Gamma) frailty models with three different sample sizes at two censoring situations. Furthermore, we compare the performance of the two proposed profile and non-profile MM algorithms with the existing estimation approach by coxph function of the survival R package under gamma and log-normal frailty models since the coxph function of the survival R package only allows gamma and log-normal frailty. The computation time of MM algorithm can be improved by using simple off-the-shelf accelerators ([31,32]); here, we implement the accelerated profile MM and non-profile MM algorithms with the squared iterative method (SqS1). We set the stopping criterion of iteration as follows.

$$\begin{array}{c}\hfill \frac{\left|\ell ({\mathit{\alpha }}^{(k+1)}|Y_{obs})-\ell ({\mathit{\alpha }}^{(k)}|Y_{obs})\right|}{\left|\ell ({\mathit{\alpha }}^{(k)}|Y_{obs})\right|+1}<{10}^{-6}.\end{array}$$

Based on 500 replications, the average values of estimated frailty and regression parameters (MLE), their biases (Bias), their empirical standard deviations (SD), and run times (T) based on three estimation methods are summarized in

Table 1. Simulation results for Log-normal frailty model in Example 1 with different sample sizes at 15% censoring.

Par.	Scenario 1:$(\mathit{B},\mathit{M})=(\mathbf{15},\mathbf{20})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		8.5716			5.7928			0.0586
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.3438	0.0938	0.3359	0.2925	0.0425	0.3173	0.2955	0.0455	0.3458
${\beta }_1$	−2.2157	−0.2157	0.5034	−2.0969	-0.0969	0.5175	−2.1697	−0.1697	0.5842
${\beta }_5$	−2.1316	−0.1316	0.5725	−2.1392	-0.1392	0.5567	−2.2005	−0.2005	0.5747
${\beta }_{10}$	−1.0976	−0.0976	0.5523	−1.0863	−0.0863	0.5232	−1.0854	−0.0854	0.5794
${\beta }_{15}$	1.0720	0.0720	0.5160	1.0126	0.0126	0.5185	1.1232	0.1232	0.5392
${\beta }_{20}$	2.1871	0.1871	0.5584	2.0966	0.0966	0.5286	2.0560	0.0560	0.5217
${\beta }_{25}$	3.2619	0.2619	0.5576	3.1412	0.1412	0.5455	3.2396	0.2396	0.6022
${\beta }_{30}$	3.2156	0.2156	0.5775	3.1507	0.1507	0.5468	3.2781	0.2781	0.5669
Par.	Scenario 2: $(\mathit{B},\mathit{M})=(\mathbf{30},\mathbf{13})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		5.7458			7.5817			0.0893
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.3222	0.0722	0.3006	0.2520	0.0020	0.1099	0.2813	0.0313	0.1065
${\beta }_1$	−2.1110	−0.1110	0.4399	−2.0930	−0.0930	0.4531	−2.0783	−0.0783	0.4743
${\beta }_5$	−2.1236	−0.1236	0.4634	−2.1208	−0.1208	0.4615	−2.1232	−0.1232	0.4384
${\beta }_{10}$	−1.0788	−0.0788	0.4491	−1.0528	−0.0528	0.4517	−1.0598	−0.0598	0.4877
${\beta }_{15}$	1.0507	0.0507	0.4467	0.9944	−0.0056	0.4416	1.0457	0.0457	0.4539
${\beta }_{20}$	2.0663	0.0663	0.4679	2.0421	0.0421	0.4665	2.1415	0.1415	0.4872
${\beta }_{25}$	3.1734	0.1734	0.4824	3.0586	0.0586	0.4499	3.2054	0.2054	0.4726
${\beta }_{30}$	3.1833	0.1833	0.4574	3.0666	0.0666	0.4371	3.1279	0.1279	0.4558
Par.	Scenario 3: $(\mathit{B},\mathit{M})=(\mathbf{50},\mathbf{10})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		6.2193			11.4968			0.1281
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.2683	0.0183	0.2427	0.2478	−0.0022	0.0862	0.2866	0.0366	0.0884
${\beta }_1$	−2.0630	−0.0630	0.3870	−2.0491	−0.0491	0.3831	−2.0548	−0.0548	0.4021
${\beta }_5$	−2.0668	−0.0668	0.3975	−2.0430	−0.0430	0.3903	−2.0340	−0.0340	0.3716
${\beta }_{10}$	−1.0286	−0.0286	0.3875	−1.0327	−0.0327	0.3799	−1.0583	−0.0583	0.4058
${\beta }_{15}$	1.0613	0.0613	0.4001	1.0131	0.0131	0.3915	1.0136	0.0136	0.3867
${\beta }_{20}$	2.0946	0.0946	0.3984	2.0315	0.0315	0.3853	2.0697	0.0697	0.4497
${\beta }_{25}$	3.1535	0.1535	0.4104	3.0739	0.0739	0.3973	3.1178	0.1178	0.4147
${\beta }_{30}$	3.1266	0.1266	0.4187	3.0378	0.0378	0.4088	3.1216	0.1216	0.3875

,

Table 2. Simulation results for Log-normal frailty model in Example 1 with different sample sizes at 30% censoring.

Par.	Scenario 1:$(\mathit{B},\mathit{M})=(\mathbf{15},\mathbf{20})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		4.8084			5.7248			0.0541
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.2712	0.0212	0.1446	0.2491	−0.0009	0.1434	0.3022	0.0522	0.1507
${\beta }_1$	−2.2181	−0.2181	0.6137	−2.1211	−0.1211	0.5712	−2.1928	−0.1928	0.6084
${\beta }_5$	−2.1717	−0.1717	0.6270	−2.1431	−0.1431	0.6230	−2.2133	−0.2133	0.5920
${\beta }_{10}$	−1.1136	−0.1136	0.6001	−1.0786	−0.0786	0.5881	−1.0936	−0.0936	0.6158
${\beta }_{15}$	1.0926	0.0926	0.6117	1.0452	0.0452	0.5603	1.0499	0.0499	0.6049
${\beta }_{20}$	2.2211	0.2211	0.6234	2.0423	0.0423	0.5667	2.2068	0.2068	0.6567
${\beta }_{25}$	3.2667	0.2667	0.6156	3.1369	0.1369	0.6192	3.2481	0.2481	0.5558
${\beta }_{30}$	3.2490	0.2490	0.6192	3.1452	0.1452	0.6355	3.2599	0.2599	0.5502
Par.	Scenario 2: $(\mathit{B},\mathit{M})=(\mathbf{30},\mathbf{13})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		5.56726			7.7258			0.0818
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.3189	0.0689	0.9271	0.2496	−0.0004	0.1150	0.2749	0.0249	0.1077
${\beta }_1$	−2.0826	−0.0826	0.4974	−2.0987	−0.0987	0.4898	−2.0761	−0.0761	0.5023
${\beta }_5$	−2.1135	−0.1135	0.5391	−2.0965	−0.0965	0.4950	−2.0974	−0.0974	0.5200
${\beta }_{10}$	−1.0764	−0.0764	0.5065	−1.0743	−0.0743	0.5180	−1.0510	−0.0510	0.5055
${\beta }_{15}$	1.0654	0.0654	0.5137	1.0433	0.0433	0.5028	1.0739	0.0739	0.5766
${\beta }_{20}$	2.1250	0.1250	0.5033	2.0525	0.0525	0.4933	2.1029	0.1029	0.5137
${\beta }_{25}$	3.1731	0.1731	0.5236	3.0821	0.0821	0.5105	3.1316	0.1316	0.5783
${\beta }_{30}$	3.1966	0.1966	0.5207	3.0716	0.0716	0.5034	3.2236	0.2236	0.5268
Par.	Scenario 3: $(\mathit{B},\mathit{M})=(\mathbf{50},\mathbf{10})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		4.7591			8.4607			0.1188
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.2729	0.0229	0.4236	0.2434	−0.0066	0.0865	0.2792	0.0292	0.0966
${\beta }_1$	−2.0917	−0.0917	0.3814	−2.0714	−0.0714	0.3741	−2.1494	−0.1494	0.4376
${\beta }_5$	−2.0692	−0.0692	0.4076	−2.0495	−0.0495	0.3990	−2.0617	−0.0617	0.4617
${\beta }_{10}$	−1.0320	−0.0320	0.3822	−1.0341	−0.0341	0.3725	−1.0714	−0.0714	0.4174
${\beta }_{15}$	1.0606	0.0606	0.3928	1.0163	0.0163	0.3843	1.0991	0.0991	0.3947
${\beta }_{20}$	2.0530	0.0530	0.3851	1.9865	−0.0135	0.3773	2.0780	0.0780	0.4779
${\beta }_{25}$	3.1058	0.1058	0.4000	3.0158	0.0158	0.3855	3.1531	0.1531	0.4395
${\beta }_{30}$	3.1071	0.1071	0.3849	3.0182	0.0182	0.3752	3.1250	0.1250	0.4593

,

Table 3. Simulation results for Gamma frailty model in Example 1 with different sample sizes at 15% censoring.

Par.	Scenario 1:$(\mathit{B},\mathit{M})=(\mathbf{15},\mathbf{20})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		22.3576			18.1613			0.0846
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	2.0584	0.0584	0.6494	2.1028	0.1028	0.7066	2.1197	0.1197	0.7023
${\beta }_1$	−2.1343	−0.1343	0.5452	−2.1155	−0.1155	0.4974	−2.1027	−0.1027	0.6415
${\beta }_5$	−2.1700	−0.1700	0.5426	−2.1303	−0.1303	0.4923	−2.1633	−0.1633	0.5354
${\beta }_{10}$	−1.0561	−0.0561	0.5395	−1.0550	−0.0550	0.5505	−1.0768	−0.0768	0.6086
${\beta }_{15}$	1.0515	0.0515	0.5749	1.0073	0.0073	0.5046	1.0883	0.0883	0.6006
${\beta }_{20}$	2.1385	0.1385	0.542	2.1311	0.1311	0.5639	2.1685	0.1685	0.5608
${\beta }_{25}$	3.1702	0.1702	0.5815	3.1928	0.1928	0.5893	3.2140	0.2140	0.5916
${\beta }_{30}$	3.2185	0.2185	0.5623	3.1268	0.1268	0.6005	3.2433	0.2433	0.6320
Par.	Scenario 2: $(\mathit{B},\mathit{M})=(\mathbf{30},\mathbf{13})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		23.2786			41.7191			0.1438
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	2.0934	0.0934	0.4866	2.0270	0.0270	0.4817	2.151	0.1510	0.6232
${\beta }_1$	−2.1266	−0.1266	0.4648	−2.1027	−0.1027	0.4600	−2.1189	−0.1189	0.5060
${\beta }_5$	−2.1197	−0.1197	0.4763	−2.0886	−0.0886	0.4549	−2.1842	−0.1842	0.4881
${\beta }_{10}$	−1.0328	−0.0328	0.4703	−1.0685	−0.0685	0.4789	−1.0671	−0.0671	0.4818
${\beta }_{15}$	1.0604	0.0604	0.4871	1.0212	0.0212	0.4847	1.0956	0.0956	0.4446
${\beta }_{20}$	2.1313	0.1313	0.4696	2.0185	0.0185	0.4745	2.0961	0.0961	0.4848
${\beta }_{25}$	3.1788	0.1788	0.4695	3.0747	0.0747	0.4954	3.1749	0.1749	0.4888
${\beta }_{30}$	3.1934	0.1934	0.5120	3.1287	0.1287	0.4953	3.1908	0.1908	0.5070
Par.	Scenario 3: $(\mathit{B},\mathit{M})=(\mathbf{50},\mathbf{10})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		27.4950			47.8408			0.2415
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	2.0792	0.0792	0.4091	2.0268	0.0268	0.3890	2.0905	0.0905	0.4287
${\beta }_1$	−2.0896	−0.0896	0.4017	−2.0466	−0.0466	0.4406	−2.1372	−0.1372	0.4366
${\beta }_5$	−2.0619	−0.0619	0.4060	−2.0388	−0.0388	0.4199	−2.1541	−0.1541	0.463
${\beta }_{10}$	−1.0431	−0.0431	0.4186	−1.0190	−0.0190	0.4047	−1.0672	−0.0672	0.4016
${\beta }_{15}$	1.0739	0.0739	0.4190	1.0363	0.0363	0.4025	1.0006	0.0006	0.4269
${\beta }_{20}$	2.0655	0.0655	0.4141	2.0694	0.0694	0.4272	2.0781	0.0781	0.4213
${\beta }_{25}$	3.1542	0.1542	0.3919	3.0680	0.0680	0.4461	3.1171	0.1171	0.4331
${\beta }_{30}$	3.1085	0.1085	0.4201	3.0707	0.0707	0.4800	3.1198	0.1198	0.4200

,

Table 4. Simulation results for Gamma frailty model in Example 1 with different sample sizes at 30% censoring.

Par.	Scenario 1:$(\mathit{B},\mathit{M})=(15,20)$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		18.1609			27.0720			0.0601
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	2.0325	0.0325	0.8620	2.0844	0.0844	0.7522	2.0332	0.0332	0.6807
${\beta }_1$	−1.9628	0.0372	0.8105	−2.2108	−0.2108	0.6094	−2.1984	−0.1984	0.5954
${\beta }_5$	−1.9076	0.0924	0.8051	−2.1330	−0.1330	0.5744	−2.1072	−0.1072	0.6045
${\beta }_{10}$	−1.0357	−0.0357	0.5701	−1.0732	−0.0732	0.5423	−1.1172	−0.1172	0.5989
${\beta }_{15}$	0.9478	−0.0522	0.6464	1.0731	0.0731	0.5596	1.0585	0.0585	0.5717
${\beta }_{20}$	1.9482	−0.0518	0.7783	2.0924	0.0924	0.5450	2.1451	0.1451	0.6499
${\beta }_{25}$	2.9937	−0.0063	1.0379	3.1667	0.1667	0.6197	3.2788	0.2788	0.5763
${\beta }_{30}$	2.9867	−0.0133	1.0485	3.2437	0.2437	0.6203	3.2631	0.2631	0.5981
Par.	Scenario 2: $(\mathit{B},\mathit{M})=(\mathbf{30},\mathbf{13})$
	Profile MM approach			Non-Profile MM Approach			Coxph
T		16.7308			30.7202			0.0995
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	2.0880	0.0880	0.5207	2.0962	0.0962	0.4792	2.1913	0.1913	0.6610
${\beta }_1$	−2.1346	−0.1346	0.5360	−2.1102	−0.1102	0.4817	−2.1127	−0.1127	0.5501
${\beta }_5$	−2.0743	−0.0743	0.5299	−2.1301	−0.1301	0.5327	−2.0796	−0.0796	0.5362
${\beta }_{10}$	−1.0846	−0.0846	0.5289	−1.0862	−0.0862	0.5136	−1.0690	−0.0690	0.4962
${\beta }_{15}$	1.0980	0.0980	0.5242	1.0001	0.0001	0.5169	1.0932	0.0932	0.5306
${\beta }_{20}$	2.1329	0.1329	0.5216	2.0354	0.0354	0.5278	2.1117	0.1117	0.5432
${\beta }_{25}$	3.2089	0.2089	0.5183	3.1365	0.1365	0.5155	3.2383	0.2383	0.5587
${\beta }_{30}$	3.1702	0.1702	0.5345	3.1695	0.1695	0.5495	3.2301	0.2301	0.5380
Par.	Scenario 3: $(\mathit{B},\mathit{M})=(\mathbf{50},\mathbf{10})$
	Profile MM Approach			Non-Profile MM Approach			Coxph
T		16.3097			29.1778			0.1647
	MLE	Bias	SD	MLE	Bias	SD	MLE	Bias	SD
$\theta$	2.0999	0.0999	0.4480	2.0462	0.0462	0.3831	2.0948	0.0948	0.3700
${\beta }_1$	−2.1151	−0.1151	0.4605	−2.0508	−0.0508	0.4187	−2.1324	−0.1324	0.4554
${\beta }_5$	−2.0767	−0.0767	0.4197	−2.0876	−0.0876	0.4615	−2.1001	−0.1001	0.4394
${\beta }_{10}$	−1.0963	−0.0963	0.4629	−1.0712	−0.0712	0.4393	−1.0513	−0.0513	0.4306
${\beta }_{15}$	1.0443	0.0443	0.4313	0.9910	−0.0090	0.4234	2.0347	0.0347	0.4476
${\beta }_{20}$	2.0813	0.0813	0.4775	2.0330	0.0330	0.4630	2.1039	0.1039	0.4538
${\beta }_{25}$	3.1089	0.1089	0.4848	3.0978	0.0978	0.4660	3.1829	0.1829	0.4484
${\beta }_{30}$	3.1216	0.1216	0.4805	3.1462	0.1462	0.5005	3.1834	0.1834	0.5308

,

Table 5. Simulation results for Inverse Gaussian model in Example 1 with different sample sizes at 15% censoring.

Par.	Scenario 1:$(\mathit{B},\mathit{M})=(\mathbf{15},\mathbf{20})$
	Profile MM Approach			Non-Profile MM Approach
T		3.0369			7.9534
	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.7060	−0.2940	0.6920	0.7130	−0.2870	0.6511
${\beta }_1$	−2.1516	−0.1516	0.5664	−2.1002	−0.1002	0.4646
${\beta }_5$	−2.1783	−0.1783	0.5366	−2.1062	−0.1062	0.5173
${\beta }_{10}$	−1.1028	−0.1028	0.5526	−1.0881	−0.0881	0.4719
${\beta }_{15}$	1.1284	0.1284	0.5523	0.9661	−0.0339	0.5069
${\beta }_{20}$	2.2162	0.2162	0.5225	1.9742	−0.0258	0.4936
${\beta }_{25}$	3.2591	0.2591	0.5269	2.9877	−0.0123	0.5332
${\beta }_{30}$	3.2139	0.2139	0.5553	3.0247	0.0247	0.5250
Par.	Scenario 2:$(\mathit{B},\mathit{M})=(\mathbf{30},\mathbf{13})$
	Profile MM Approach			Non-Profile MM Approach
T		5.2872			27.2494
	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.7537	−0.2463	0.5623	0.9178	−0.0822	0.4438
${\beta }_1$	−2.1515	−0.1515	0.4592	−2.1441	−0.1441	0.4455
${\beta }_5$	−2.1300	−0.1300	0.4646	−2.1146	−0.1146	0.4907
${\beta }_{10}$	−1.0873	−0.0873	0.4743	−1.0830	−0.0830	0.4316
${\beta }_{15}$	1.0681	0.0681	0.4755	1.0223	0.0223	0.4445
${\beta }_{20}$	2.1454	0.1454	0.4600	2.0591	0.0591	0.4625
${\beta }_{25}$	3.2323	0.2323	0.4972	3.1031	0.1031	0.4807
${\beta }_{30}$	3.2223	0.2223	0.4698	3.0855	0.0855	0.4907
Par.	Scenario 3:$(\mathit{B},\mathit{M})=(\mathbf{50},\mathbf{10})$
	Profile MM Approach			Non-Profile MM Approach
T		54.9905			36.7796
	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.8461	−0.1539	0.4143	0.9724	−0.0276	0.3678
${\beta }_1$	−2.0933	−0.0933	0.4056	−2.0688	−0.0688	0.4135
${\beta }_5$	−2.1509	−0.1509	0.4375	−2.0664	−0.0664	0.4318
${\beta }_{10}$	−1.0878	−0.0878	0.3992	−1.0790	−0.0790	0.3824
${\beta }_{15}$	1.0437	0.0437	0.4217	0.9728	−0.0272	0.4111
${\beta }_{20}$	2.1137	0.1137	0.4077	2.0257	0.0257	0.4005
${\beta }_{25}$	3.1511	0.1511	0.4311	3.0709	0.0709	0.4129
${\beta }_{30}$	3.1846	0.1846	0.4137	3.0540	0.0540	0.4219

and

Table 6. Simulation results for Inverse Gaussian model in Example 1 with different sample sizes at 30% censoring.

Par.	Scenario 1:$(\mathit{B},\mathit{M})=(\mathbf{15},\mathbf{20})$
	Profile MM Approach			Non-Profile MM Approach
T		9.7444			16.1049
	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.5343	−0.4657	0.6665	0.7017	−0.2983	0.6431
${\beta }_1$	−2.1693	−0.1693	0.6153	−2.1795	−0.1795	0.6054
${\beta }_5$	−2.2288	−0.2288	0.5688	−2.1677	−0.1677	0.6601
${\beta }_{10}$	−1.0228	−0.0228	0.5690	−1.1006	−0.1006	0.5843
${\beta }_{15}$	1.1635	0.1635	0.5542	1.0577	0.0577	0.6157
${\beta }_{20}$	2.2779	0.2779	0.6479	2.0902	0.0902	0.6168
${\beta }_{25}$	3.2790	0.2790	0.6760	3.1711	0.1711	0.6672
${\beta }_{30}$	3.3969	0.3969	0.6867	3.1939	0.1939	0.6412
Par.	Scenario 2: $(\mathit{B},\mathit{M})=(\mathbf{30},\mathbf{13})$
	Profile MM Approach			Non-Profile MM Approach
T		22.7787			22.0124
	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.6564	−0.3436	0.5263	1.0203	0.0203	0.5297
${\beta }_1$	−2.1801	−0.1801	0.5395	−2.1915	−0.1915	0.4992
${\beta }_5$	−2.0213	−0.0213	0.4727	−2.1508	−0.1508	0.5142
${\beta }_{10}$	−1.1390	−0.1390	0.5202	−1.0678	−0.0678	0.5140
${\beta }_{15}$	1.0824	0.0824	0.5276	0.9845	−0.0155	0.5005
${\beta }_{20}$	2.1354	0.1354	0.4639	2.0252	0.0252	0.5088
${\beta }_{25}$	3.1809	0.1809	0.6421	3.1362	0.1362	0.5233
${\beta }_{30}$	3.2375	0.2375	0.5542	3.0828	0.0828	0.5183
Par.	Scenario 3: $(\mathit{B},\mathit{M})=(\mathbf{50},\mathbf{10})$
	Profile MM Approach			Non-Profile MM Approach
T		20.9400			32.6620
	MLE	Bias	SD	MLE	Bias	SD
$\theta$	0.8502	−0.1498	0.3856	0.9809	−0.0191	0.3556
${\beta }_1$	−2.1046	−0.1046	0.4252	−2.0688	−0.0688	0.4124
${\beta }_5$	−2.1536	−0.1536	0.4224	−2.1123	−0.1123	0.4065
${\beta }_{10}$	−1.0597	−0.0597	0.3893	−1.0542	−0.0542	0.3806
${\beta }_{15}$	1.0692	0.0692	0.3973	1.0067	0.0067	0.3841
${\beta }_{20}$	2.1129	0.1129	0.4162	2.0184	0.0184	0.4005
${\beta }_{25}$	3.1770	0.1770	0.4159	3.0489	0.0489	0.4005
${\beta }_{30}$	3.2102	0.2102	0.4277	3.0832	0.0832	0.4120

. In general, with a sample size increase, the biases and empirical standard deviations of almost all parameters become smaller. Small sample size $(B,M)=(15,20)$ always causes obvious biases for frailty parameters in Log-normal and Inverse Gaussian frailty models because the number of clusters is too small. Moreover, a larger censoring proportion usually results in greater empirical standard deviations for almostly all cases. From Table 1, Table 2, Table 3 and Table 4, we also observe that the existing estimation approach using coxph function is the fastest among the three methods since the survival R package is optimised using the C language. In terms of estimation accuracy, the non-profile MM algorithm performs the best for both frailty and regression parameters, almost always exhibiting smallest biases and empirical standard deviations in all situations. Even in a small sample size such as $(B,M)=(15,20)$, the non-profile MM estimates for frailty parameter still perform well, especially for Log-normal frailty and Gamma frailty models.

Example 2.

In this example, we simulated 200 realizations consisting of $B=50$ clusters and $M=6$ subjects in each cluster from the frailty model:

$$\begin{array}{c}\hfill \lambda (t|{\mathbf{X}}_{ij},{\omega }_i)={\omega }_i{\lambda }_0(t)exp\left\{{\mathbf{X}}_{ij}^{\top }\mathit{\beta }\right\},\end{array}$$

(28)

where the frailty terms are ${\omega }_i\sim Gamma(1/\theta ,1/\theta )$ with $\theta =0.5$ for $i=1,\dots ,B$, $\mathit{\beta }={(1,3,0_{46}^{\top },2,4)}^{\top },q=50,{\lambda }_0(t)=5$. The $X_i$ were marginally standard normal, and the correlation between $X_i$ and $X_j$ is ${\varrho }^{|i-j|}$ with $\varrho =0.25,0.75$, respectively. The censoring times were generated from uniform distribution to yield a censoring proportion around $30\%$. In this simulation experiment, the utility of our proposed profile MM and non-profile MM estimation methods for regularized estimation was illustrated in a sparse high-dimensional regression model (21) with three different penalties ($L_1$, MCP, and SCAD). The model error (ME) and relative model error (RME), which is the ratio of the model error of the regularized estimator and that of the ordinary maximum likelihood estimator, are calculated to evaluate the estimation accuracy. Based on 200 replications, we report the median of relative model errors (MRME) and the average number of correctly and incorrectly identified zero coefficients in

Table 7. The median of relative model errors for gamma frailty model by $L_1$, MCP, and SCAD penalties with sample size $(B,M)=(50,6)$ based on 200 realizations in Example 2.

	MRME	Zeros		MRME	Zeros
		Correct	Incorrect		Correct	Incorrect
Penalty	$\mathit{\varrho }=\mathbf{0.25}$			$\mathit{\varrho }=\mathbf{0.75}$
Profile MM method
$L_1$	0.159	45.775	0	0.163	44.39	0
MCP ($\gamma =3$)	0.091	46	0	0.066	46	0
SCAD ($\gamma =3.7$)	0.143	45.915	0	0.101	45.965	0
Non-profile MM method
$L_1$	0.089	45.74	0	0.188	44.225	0
MCP ($\gamma =3$)	0.051	46	0	0.085	46	0
SCAD ($\gamma =3.7$)	0.077	45.885	0	0.106	45.88	0

in which the column labeled “Correct" presents the average restricted only to the true zero coefficients, while the column labeled “Incorrect" depicts the average of coefficients erroneously set to 0. From Table 4, we can observe that the proposed profile and non-profile MM algorithms mesh very well with MCP and SCAD penalties and provide good results in both parameter estimation and variable selection, especially for MCP penalty. For $L_1$ penalty, we observe that it meshes well with profile and non-profile MM algorithms under lower correlation case at $\varrho =0.25$ but tends to yield biased estimates and inaccurate variable selection results when the correlation between $X_i$ and $X_j$ becomes stronger. Furthermore, we report the average values of estimated parameters (MLE), their biases (Bias), and their empirical standard deviations (SD) under MCP and SCAD penalties based on 200 repetition in

Table 8. The average values of estimated parameters (MLE), their biases (Bias), and their empirical standard deviations (SD) under MCP and SCAD penalties in Example 2.

Penalty	Par.	Profile MM Approach			Non-Profile MM Approach
		MLE	Bias	SD	MLE	Bias	SD
MCP (ϱ = 0.25)	$\theta$	0.4835	$-0.0165$	0.1120	0.4820	$-0.0180$	0.1398
	${\beta }_1$	$0.9993$	$-0.0007$	0.0950	0.9988	$-0.0012$	0.1028
	${\beta }_2$	$3.0090$	$0.0090$	0.1984	2.9856	$-0.0144$	0.1971
	${\beta }_{49}$	$1.9940$	$-0.0060$	0.1513	1.9791	$-0.0209$	0.1512
	${\beta }_{50}$	$3.9697$	$-0.0303$	0.2498	3.9868	$-0.0132$	0.2566
SCAD (ϱ = 0.25)	$\theta$	0.4967	$-0.0033$	0.1267	0.4800	$-0.0200$	0.1335
	${\beta }_1$	1.0133	0.0133	0.0960	0.9864	$-0.0136$	0.0975
	${\beta }_2$	3.0167	0.0167	0.2023	2.9467	$-0.0532$	0.1644
	${\beta }_{49}$	2.0183	0.0183	0.1417	1.9816	$-0.0184$	0.1255
	${\beta }_{50}$	4.0315	0.0315	0.2657	3.9300	$-0.0700$	0.2334
MCP (ϱ = 0.75)	$\theta$	0.4886	$-0.0114$	0.1258	0.4842	$-0.0158$	0.1480
	${\beta }_1$	0.9774	$-0.0226$	0.1331	1.0071	0.0071	0.1369
	${\beta }_2$	3.0252	0.0252	0.2020	2.9743	$-0.0257$	0.2570
	${\beta }_{49}$	2.0083	0.0083	0.1831	1.9857	$-0.0143$	0.2203
	${\beta }_{50}$	4.0134	0.0134	0.2606	4.0026	0.0026	0.3055
SCAD (ϱ = 0.75)	$\theta$	0.4869	$-0.0131$	0.1277	0.4839	$-0.0161$	0.1341
	${\beta }_1$	1.0119	0.0119	0.1364	1.0152	0.0152	0.1560
	${\beta }_2$	2.9815	$-0.0185$	0.2063	3.0035	0.0035	0.2422
	${\beta }_{49}$	2.0012	0.0012	0.2238	1.9984	$-0.0016$	0.2052
	${\beta }_{50}$	4.0119	0.0119	0.2633	4.0021	0.0021	0.2855

. It can be observed that both profile and non-profile MM methods perform similarly well in varying ϱ and different penalties, always showing small biases and empirical standard deviations.

6. Real Data Analysis

Alzheimer Disease’s (AD) is the most common dementia that causes progressive memory and other body function losses. According to Alzheimer Association, Alzheimer’s is the sixth leading cause of death in the United States. Many researches have studied the development of AD from Mild Cognitive Impairment (MCI). The authors of [33] analyzed the AUC score of whether MCI will transfer to AD using LASSO penalized logistic regression. The authors of [34] applied cognitive scores, ApoE genotype, and CSF biomarkers to predict the transition time from MCI to AD. Moreover, there are research studies such as [35] that used image data to analyze the transition time to dementia. In this paper, we will apply clinical data and selected SNP genotypes to conduct feature selection under three different frailty models (Gamma, Inverse Gaussian, and Log-normal) using SCAD and MCP penalties. The dataset was obtained from ADNI database (adni.loni.usc.edu). In this dataset, 267 people were recorded as cognitively normal (CN) during the first visit. Among these people, 78 of them were diagnosed with MCI before the last visit. Eventually, we observed 22 people transferred from the MCI stage to dementia. To predict the these two conversion times, 19 clinical predictors were applied mainly based on the age, marriage status, education level, and test score wuch as Alzheimer’s Disease Assessment Scale cognitive subscale (ADAS-Cog) and Functional Activities Questionnaire (FAQ). Other than these predictors, we also included genotypes of SNP from GWAS such as ApoE-$\epsilon 4$, which have a relationship with early onset dementia or late-onset dementia; 132 covariates are selected for model training.

In general, using the same notation from simulation, we have $(B,M)=(276,2)$. Individuals are independent, and two events of the same individual are grouped into the same cluster that share the same frailty. The dataset contains the following information. For individual ($i=1,2,\cdots ,276$) and event ($j=1,2$), $t_{ij}$ is minimum value of event time $Y_{ij}$ and censoring time $C_{ij}$. As shown by

Table 9. Illustration of survival years and censoring date for some patients from the ADNI dataset.

Patient ID	From CN to MCI
	First Observation Date of State CN	Date of Transition to MCI	Last Date of Observation	${\mathit{t}}_{\mathit{ij}}$ (yr)	${\mathit{\delta }}_{\mathit{ij}}$
011_S_0002	8 September 2005	26 September 2012	18 October 2017	7.05	1
011_S_0021	24 October 2005	-	27 November 2017	12.10	0
100_S_0035	8 November 2005	8 December 2010	8 December 2010	5.08	1
131_S_0123	7 February 2006	23 February 2012	10 February 2016	6.05	1
127_S_0259	28 March 2006	9 April 2014	22 September 2017	8.04	1
Patient ID	From MCI to Dementia
	First Observation Date of State MCI	Date of Transition to Dementia	Last Date of Observation	${\mathit{t}}_{\mathit{ij}}$ (yr)	${\mathit{\delta }}_{\mathit{ij}}$
011_S_0002	26 September 2012	-	18 October 2017	5.06	0
011_S_0021	-	-	-	-	-
100_S_0035	8 December 2010	-	8 December 2010	0.00	0
131_S_0123	23 February 2012	12 February 2014	10 February 2016	1.97	1
127_S_0259	9 April 2014	10 April 2015	22 September 2017	1.00	1

, censoring time $C_{ij}$ in this dataset is the time difference between the first observation date of certain event (CN or MCI) and last observation date of this patient. The event time is the state transition time, and $Y_{i1}$ is the time difference between the first observation date of this patient in state CN and date of transition from CN to MCI. If no transition is observed, censoring occurs where ${\delta }_{i1}=0$ and $t_{i1}=C_{i1}$, otherwise, ${\delta }_{i1}=1$ and $t_{i1}=Y_{i1}$. Similarly, $Y_{i2}$ is the time difference between the first observation date of this patient in state MCI and the date of transition from MCI to Dementia. If no transition is observed, ${\delta }_{i2}=0$ and $t_{i2}=C_{i2}$, otherwise, ${\delta }_{i2}=1$ and $t_{i2}=Y_{i2}$. ${\mathbf{X}}_{ij}^{\top }=(X_{ij1},\dots ,X_{ijp})$ denote the potential covariates (such as SNP genotypes and clinical predictors) where $p=132$. According to data description, it is reasonable to assume that censoring time $C_{ij}$ is independent of event time $Y_{ij}$. Moreover, the two events from each individual are considered to be highly associated due to common genetic factors and/or certain habits. Thus, the frailty model is suitable for analyzing this dataset.

Three frailty models are applied for variable selection using SCAD and MCP penalties. According to the results presented by

Table 10. Training results for Gamma Frailty Model.

	Profile MM		Non-Profile MM
Penalty	BIC	No. of Non-Zero $\mathit{\beta }$	BIC	No. of Non-Zero $\mathit{\beta }$
MCP	1154.70	70	1154.75	70
SCAD	1150.59	69	1150.82	69

,

Table 11. Training results for Inverse Gaussian Model.

	Profile MM		Non-Profile MM
Penalty	BIC	No. of Non-Zero $\mathit{\beta }$	BIC	No. of Non-Zero $\mathit{\beta }$
MCP	1198.76	27	1198.82	27
SCAD	1203.97	27	1204.50	27

and

Table 12. Training results for Log-normal Frailty Model.

	Profile MM		Non-Profile MM
Penalty	BIC	No. of Non-Zero $\mathit{\beta }$	BIC	No. of Non-Zero $\mathit{\beta }$
MCP	1149.10	71	1148.86	72
SCAD	1149.10	71	1148.86	71

, around 70 covariates were selected for Gamma Frailty and Log-normal Frailty model, and 27 covariates were selected by an Inverse Gaussian model. Similar covairates were selected by SCAD or MCP penalty for both Profile MM and Non-profile MM under the same model, which verifies the result from simulation study. For selected covariates, the ApoE-$\epsilon 4$ is selected as a negative effect to the survival time as expected. In addition to ApoE-$\epsilon 4$, single SNPs such as rs2333227 and rs669 from the MPO and A2M gene are associated with the development of Alzheimer’s Disease. We can observe similar results from [36]. However, many selected SNPs are not recognized from other studies, and it needs to be analyzed further whether these SNPs can help to identify the risk of Alzheimer’s.

We use the model with the lowest BIC score for prediction by constructing prognostic index in order to see whether the selected covariates can help identify the individual with low risk or high risk. Let $\mathbf{X}$ be a collection of ${\mathbf{X}}_{ij}^{\top }$. The full dataset is divided into 10 parts (${\mathbf{X}}^{(k)}$$(k=1,2,\dots ,10)$). Similarly to the cross-validation process, without the part of dataset ${\mathbf{X}}^{(k)}$, covariates are selected using BIC criteria and the vector of parameters ${\widehat{\beta }}_k$ is estimated correspondingly. The risk factor ${\widehat{r}}_k={\mathbf{X}}^{(k)}{\widehat{\beta }}_k$$(k=1,2,\dots ,10)$ is calculated using dataset ${\mathbf{X}}^{(k)}$ without considering the individuals’ Latent variable. Individuals with estimated risk factors below median are classified into to low risk group, and and individual with estimated risk above median is classified into the high risk group. The survival plots for estimated high risk group and low risk group are constructed based on their true survival times. Figure 1 shows that the selected model has conducted a good estimation for the transition time from stage CN to stage MCI. Due to the lack of data for the second transition, their confidence intervals are overlapped. Therefore, the prediction result, especially from stage CN to stage MCI, has shown good estimation using the corresponding model, and two risk groups are distinguished by predicted risk levels.

7. Discussion

The profile and non-profile MM algorithms are proposed for high-dimension regression analysis with clustered failure time data where a general frailty is used to model within-cluster dependence and the penalty such as the SCAD and MCP is used for inducing sparsity. The proposed methods can separate the high-dimensional minorizing function into a sum of univariate functions after a sequence of minorization steps. These approaches avoid matrix inversion and provide a toolkit for developing more efficient algorithms in a broad range of in statistical optimization problems. Meshing well with sparsity penalties such as SCAD or MCP, the two regularized MM algorithms are further shown to exhibit certain numerical advantages in sparse high-dimensional regression analyses. The shared frailty model only represents a special and relatively simple model among the widely used frailty models for multivariate survival data. For example, the standard shared frailty model assumes that all subjects in the same cluster share a common frailty. This assumption can be relaxed to the correlated frailty terms among subjects in the same cluster. Correlated frailty models present the limitation that shared frailty models may only be used to fit positively correlated event times. Furthermore, frailty is assumed to be time-constant. However, unobserved heterogeneity may also be time dependent, which can be explained by an unobserved random process that unfolds over time. Based on this idea, several approaches have been proposed such as diffusion processes modeling or Levy processes modeling for frailty. As an approach based on birth-death Poisson or simpler, piecewise constant, frailty models have recently been proposed. It would be worthwhile to extend the proposed two MM algorithms in these applications. Lastly, spatially correlated survival data present another important and useful setting where the proposed MM algorithms can be further extended to accommodate.