Sparse Boosting Based Machine Learning Methods for High-Dimensional Data

Abstract

In high-dimensional data, penalized regression is often used for variable selection and parameter estimation. However, these methods typically require time-consuming cross-validation methods to select tuning parameters and retain more false positives under high dimensionality. This chapter discusses sparse boosting based machine learning methods in the following high-dimensional problems. First, a sparse boosting method to select important biomarkers is studied for the right censored survival data with high-dimensional biomarkers. Then, a two-step sparse boosting method to carry out the variable selection and the model-based prediction is studied for the high-dimensional longitudinal observations measured repeatedly over time. Finally, a multi-step sparse boosting method to identify patient subgroups that exhibit different treatment effects is studied for the high-dimensional dense longitudinal observations. This chapter intends to solve the problem of how to improve the accuracy and calculation speed of variable selection and parameter estimation in high-dimensional data. It aims to expand the application scope of sparse boosting and develop new methods of high-dimensional survival analysis, longitudinal data analysis, and subgroup analysis, which has great application prospects.

1. Introduction

High-dimensional model has become very popular in statistical literature and many new machine learning techniques have been developed to deal with data with very large number of features. In the past decades, researchers have done a great deal of high-dimensional data analysis where the sample size $n$ is relatively small but the number of features $p$ under consideration is extremely large. It is widely known that including irrelevant predictors in the statistical model may result in unstable estimation and dreadful computing issues. Thus, variable selection is crucial to address the challenges. Among all developments, regularization procedures such as LASSO [1], smoothly clipped absolute deviation (SCAD) [2], MCP [3] and their various extensions [4, 5, 6] have been thoroughly studied and widely used to perform variable selection and estimation simultaneously in order to improve the prediction accuracy and interpretability of the statistical model. However, those penalized estimation approaches all have some tuning parameters required to be selected by computationally expensive methods like cross-validation.

In recent years, machine learning methods such as boosting have become very prominent for high-dimensional data settings since they can improve the selection accuracy substantially and reduce the chance of including irrelevant features. The original boosting algorithms were proposed by Schapire [7] which is an ensemble method that iteratively combines weaker learners to minimize the expected loss. The major difference among different boosting algorithms is the loss function. For example, AdaBoost [8] has the exponential loss function, L2 boosting [9] has the squared error loss function, sparse boosting [10] has the penalized loss function and HingeBoost [11] has the weighted hinge loss function. Recently, more various versions of boosting algorithms have been proposed. See, for example, Bühlmann and Hothorn [12] for the twin boosting; Komori and Eguchi [13] for the pAUCBoost; Wang [14] for the twin HingeBoost; Zhao [15] for the GSBoosting and Yang and Zou [16] for the ER-Boost. Besides these extensions, much effort has been made in understanding the advantages of boosting such as relatively lower over-fitting risk, smaller computational cost, and simpler adjustment to include additional constraints.

In this chapter we review some sparse boosting based methods for the following high-dimensional problems based on three research papers. First, a sparse boosting method to select important biomarkers is studied for the right censored survival data with high-dimensional biomarkers [17]. Then, a two-step sparse boosting to carry out the variable selection and the model-based prediction is studied for the high-dimensional longitudinal observations measured repeated over time [18]. Finally, a multi-step sparse boosting method to identify patient subgroups that exhibit different treatment effects is studied for the high-dimensional dense longitudinal observations [19]. This chapter intends to solve the problem of how to improve the accuracy and calculation speed of variable selection and parameter estimation in high-dimensional data. It aims to expand the application scope of sparse boosting and develop new methods of high-dimensional survival analysis, longitudinal data analysis, and subgroup analysis, which has great application prospects.

The rest of the chapter is arranged as follows. In Section 2, a sparse boosting method to fit high-dimensional survival data is studied. In Section 3, a two-step sparse boosting approach to carry out variable selection and model-based prediction by fitting high-dimensional models with longitudinal data is studied. In Section 4, a subgroup identification method incorporating multi-step sparse boosting algorithm for high-dimensional dense longitudinal data is studied. Finally, Section 5 provides concluding remarks.

2. Sparse boosting for survival data

Survival time data are usually referred to time-to-event data and they are usually censored. Predicting survival time and identifying the risk factors can be very helpful for patient treatment selection, disease prevention strategy or disease management in evidence-based medicine. A well-known model in survival analysis is the Cox proportional hazards (PH) model [20] which assumes multiplicative covariate effects in the hazards function. Another popular model is the accelerated failure time (AFT) model [21] which assumes that the covariate effect is to accelerate or decelerate the life time of a disease. The coefficients in the regression model have the direct interpretation of the covariate effects on the mean survival time. Recently, researchers developed boosting methods to analyze survival data. For example, Schmid and Hothorn [22] proposed a flexible boosting method for parametric AFT models, and Wang and Wang [23] proposed Buckley-James boosting for survival data with right censoring and high dimensionality.

In this section, a sparse boosting method to fit high-dimensional varying-coefficient AFT models is presented. In particular, the sparse boosting techniques for right censored survival data is studied. In Section 2.1, the varying-coefficient AFT model for survival data is formulated and a detailed sparse boosting algorithm to fit the model is proposed. In Section 2.2, the proposed sparse boosting techniques through simulation studies is evaluated. In Section 2.3, the performance of sparse boosting via a lung cancer data example is examined.

2.1 Methodology

2.1.1 Model and estimation

Let $T_i$ and $C_i$ be the logarithm of survival time and censoring time for the $i$th subject in a random sample of size $n$ respectively. In reality $Y_i=min\left\{T_i,C_i\right\}$ and the censoring indicator ${\delta }_i=I\left(T_i\le C_i\right)$ [24] are observed. Denote ${\mathbf{X}}_i=\left(X_{i,1},\cdots ,X_{i,p-1}\right)$ to be the corresponding ($p$-$1$)-dimensional predictors such as gene expressions or biomarkers for the $i$th subject and $U_i$ to be the univariate index variable. Our observed data set $\left\{\left({\mathbf{X}}_i,{\delta }_i,Y_i,U_i\right):{\mathbf{X}}_i\in \mathrm{I}{\mathrm{R}}^{p-1},{\delta }_i\in \left\{0,1\right\},Y_i\in \mathrm{I}\mathrm{R},U_i\in \mathrm{I}\mathrm{R},i=1,2,\cdots ,n\right\}$ is an independently and identically distributed random sample from $\left(\mathbf{X},\delta ,Y,U\right)$. The varying-coefficient AFT model is:

$$T_i={\beta }_0\left(U_i\right)+\sum _{j=1}^{p-1}{X}_{i,j}{\beta }_j\left(U_i\right)+{\epsilon }_i,i=1,\dots ,n,$$E1

where ${\beta }_0\left(.\right),{\beta }_1\left(.\right),\cdots ,{\beta }_{p-1}\left(.\right)$ are the unknown varying-coefficient functions of confounder $U$ and ${\epsilon }_i$ is the random error with $E\left({\epsilon }_i|{\mathbf{X}}_i,U_i\right)=0$.

A weighted least squares estimation approach is adopted. Let $w_i$‘s be the Kaplan–Meier weights [25], which are the jumps in the Kaplan–Meier estimator computed as $w_1=\frac{{\delta }_{\left(1\right)}}{n}$ and $w_i=\frac{{\delta }_{\left(i\right)}}{n-i+1}{\prod }_{j=1}^{i-1}{\left(\frac{n-j}{n-j+1}\right)}^{{\delta }_{\left(j\right)}},i=2,\dots ,n$. Let $Y_{\left(1\right)}\le \cdots \le Y_{\left(n\right)}$ be the order statistics of $Y_{i^{\prime }}s$, ${\delta }_{\left(1\right)},\cdots ,{\delta }_{\left(n\right)}$ be the corresponding censoring indicators of the ordered $Y_{i^{\prime }}s$, and $X_{\left(1\right),j},\cdots ,X_{\left(n\right),j},j=1,\cdots ,p-1$ and $U_{\left(1\right)},\cdots ,U_{\left(n\right)}$ are defined similarly. Then the weighed least squares loss function is

$$\sum _{i=1}^n{w}_i{\left(Y_{\left(i\right)}-{\beta }_0\left(U_{\left(i\right)}\right)-\sum _{j=1}^{p-1}{X}_{\left(i\right),j}{\beta }_j\left(U_{\left(i\right)}\right)\right)}^2.$$E2

Let $B\left(.\right)={\left(B_1\left(.\right),\dots ,B_L\left(.\right)\right)}^T$ be an equal-spaced B-spline basis, where $L$ is the dimension of the basis. Under certain smoothness conditions, the Curry-Schonberg theorem [26] implies that for every smooth function ${\beta }_j\left(.\right)$, it can be approximated by

$${\beta }_j\left(.\right)\approx B^T\left(.\right){\gamma }_j,j=0,\cdots ,p-1,$$E3

where ${\gamma }_j$ is a vector of length $L$. Then the weighted least squares loss function Eq. (2) can be approximated by

$$\sum _{i=1}^n{w}_i{\left(Y_{\left(i\right)}-B^T\left(U_{\left(i\right)}\right){\gamma }_0-\sum _{j=1}^{p-1}{X}_{\left(i\right),j}{B}^T\left(U_{\left(i\right)}\right){\gamma }_j\right)}^2.$$E4

Denote by $\tilde{Y}={\left(Y_{\left(1\right)},\cdots ,Y_{\left(n\right)}\right)}^T,$$X_{\left(i\right),0}=1$ for $i=1,\cdots ,n$, ${\tilde{\mathbf{X}}}_j={\left(B\left(U_{\left(1\right)}\right)X_{\left(1\right),j},\cdots ,B\left(U_{\left(n\right)}\right)X_{\left(n\right),j}\right)}^T$, $\tilde{\mathbf{X}}=\left({\tilde{\mathbf{X}}}_0,\cdots ,{\tilde{\mathbf{X}}}_{p-1}\right)$, $W=diag\left(w_1,\cdots ,w_n\right)$ and $\gamma ={\left({\gamma }_0^T,\cdots ,{\gamma }_{p-1}^T\right)}^T$. Then the objective function Eq. (4) may be written in the following matrix form:

$${\left(\tilde{Y}-\tilde{\mathbf{X}}\gamma \right)}^{T}W\left(\tilde{Y}-\tilde{\mathbf{X}}\gamma \right).$$E5

The estimation may yield close-form solution for the coefficients when dimensionality $p$ is small or moderate. With high dimensionality the solution cannot be easily achieved. Let ${\gamma }^{\left[\widehat{K}\right]}={\left({\left({\gamma }_0^{\left[\widehat{K}\right]}\right)}^T,\cdots ,{\left({\gamma }_{p-1}^{\left[\widehat{K}\right]}\right)}^T\right)}^T$ be the estimator of $\gamma$ from sparse boosting approach with weighted square loss function Eq. (5), and $\widehat{K}$ is the estimated number of stopping iterations. Then the estimates of coefficient function are given by

$${\widehat{\beta }}_j\left(u\right)=B^T\left(u\right){\gamma }_j^{\left[\widehat{K}\right]},j=0,\cdots ,p-1.$$E6

Instead of using the regularized estimation approaches, a sparse boosting method to estimate ${\gamma }^{\left[\widehat{K}\right]}$ is presented in the following subsection.

2.1.2 Sparse boosting techniques

The key idea of sparse boosting is to replace the empirical risk function in L2 boosting with the penalized empirical risk function which is a combination of squared loss and the trace of boosting operator as a measure of boosting complexity, and then perform gradient descent in a function space iteratively. Thus sparse boosting produces sparser models compared to L2 boosting. The g-prior minimum description length (gMDL) proposed by [27] can be used as the penalized empirical risk function to estimate the update criterion in each iteration and the stopping criterion. The gMDL takes the form:

$$\begin{array}{ll}\text{gMDL}\left(\mathit{RSS},\text{trace}\left(\mathcal{B}\right)\right)& =log\left(S\right)+\frac{\text{trace}\left(\mathcal{B}\right)}{n}log\left(\frac{{\tilde{Y}}^T\tilde{Y}-\mathit{RSS}}{\text{trace}\left(\mathcal{B}\right)\times S}\right),\\ S& =\frac{\mathit{RSS}}{n-\text{trace}\left(\mathcal{B}\right)}.\end{array}$$E7

Here $\mathit{RSS}$ is the residual sum of squares and $\mathcal{B}$ is the boosting operator. The model that achieves the shortest description of data will be selected. The advantage is that it has a data-dependent penalty for each dimension since it is explicitly given as a function of data only, thus the selection of the tuning parameter can be avoided.

The sparse boosting procedure is described in details. The initial value of $\gamma$ is set to be a zero vector, i.e. ${\gamma }^{\left[k\right]}=\mathbf{0}$ for $k=0$, while in each of the $k$th iteration ($1\le k\le K$ for $K$ being the total number of iterations) only the current residual $R^{\left[k\right]}=\tilde{Y}-\tilde{\mathbf{X}}{\gamma }^{\left[k-1\right]}$ is used to regress every $j$th working element ${\tilde{\mathbf{X}}}_j,j=0,\cdots ,p-1$. The fit denoted by ${\widehat{\lambda }}_j^{\left[k\right]}$ can be obtained by minimizing the weighted squared loss function ${\left(R^{\left[k\right]}-{\tilde{\mathbf{X}}}_j\lambda \right)}^{T}W\left(R^{\left[k\right]}-{\tilde{\mathbf{X}}}_j\lambda \right)$ with respect to $\lambda$. Hence the weighted least squared estimate is ${\widehat{\lambda }}_j^{\left[k\right]}={\left[{\left({\tilde{\mathbf{X}}}_j\right)}^{T}W\left({\tilde{\mathbf{X}}}_j\right)\right]}^{-1}{\left({\tilde{\mathbf{X}}}_j\right)}^T{\mathit{WR}}^{\left[k\right]}$, the corresponding hat matrix is ${\mathcal{H}}_j=\left({\tilde{\mathbf{X}}}_j\right){\left[{\left({\tilde{\mathbf{X}}}_j\right)}^{T}W\left({\tilde{\mathbf{X}}}_j\right)\right]}^{-1}{\left({\tilde{\mathbf{X}}}_j\right)}^{T}W$ and the weighted residual sum of squares is ${\mathit{RSS}}_j^{\left[k\right]}={\left(R^{\left[k\right]}-{\tilde{\mathbf{X}}}_j{\widehat{\lambda }}_j^{\left[k\right]}\right)}^{T}W\left(R^{\left[k\right]}-{\tilde{\mathbf{X}}}_j{\widehat{\lambda }}_j^{\left[k\right]}\right)$. The selected component ${\widehat{s}}_k$ can be obtained by:

$${\widehat{s}}_k={\text{argmin}}_{0\le j\le p-1}\text{gMDL}\left({\mathit{RSS}}_j^{\left[k\right]},\text{trace}\left({\mathcal{B}}_j^{\left[k\right]}\right)\right),$$E8

where ${\mathcal{B}}_j^{\left[1\right]}={\mathcal{H}}_j$ and ${\mathcal{B}}_j^{\left[k\right]}=I-\left(I-{\mathcal{H}}_j\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k-1}}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$ for $k>1$ is the boosting operator for selecting $j$th component in the $k$th iteration. Therefore, at each iteration there is only one working component ${\tilde{\mathbf{X}}}_{{\widehat{s}}_k}$ to be chosen, and only the corresponding coefficient vector ${\gamma }_{{\widehat{s}}_k}^{\left[k\right]}$ changes, i.e. ${\gamma }_{{\widehat{s}}_k}^{\left[k\right]}={\gamma }_{{\widehat{s}}_k}^{\left[k-1\right]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_k}^{\left[k\right]}$, where $\nu$ is the step size, while all the other ${\gamma }_j^{\left[k\right]}$ for $j\ne {\widehat{s}}_k$ remain the same. This process is repeated for $K$ iterations and estimate the stopping iteration $K$ by.

$$\widehat{K}={\text{argmin}}_{1\le k\le K}\text{gMDL}\left({\mathit{RSS}}_{{\widehat{s}}_k}^{\left[k\right]},\text{trace}\left({\mathcal{B}}^{\left[k\right]}\right)\right),$$E9

where ${\mathcal{B}}^{\left[k\right]}=I-\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_k}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$.

From this sparse boosting procedure, the estimator of $\gamma$ is obtained as ${\gamma }^{\left[\widehat{K}\right]}={\left({\left({\gamma }_0^{\left[\widehat{K}\right]}\right)}^T,\cdots ,{\left({\gamma }_{p-1}^{\left[\widehat{K}\right]}\right)}^T\right)}^T$. The sparse boosting algorithm for the varying-coefficient AFT model can be summarized as follows:

Sparse Boosting Algorithm for Varying-Coefficient AFT Model.

1. Initialization. Set $k=0$ and ${\gamma }_0^{\left[k\right]}=\mathbf{0},\cdots ,{\gamma }_{p-1}^{\left[k\right]}=\mathbf{0}$ (component-wise).

2. Iteration. $k=k+1$. Compute ${\widehat{s}}_k={\text{argmin}}_{0\le j\le p-1}\text{gMDL}\left({\mathit{RSS}}_j^{\left[k\right]},\text{trace}\left({\mathcal{B}}_j^{\left[k\right]}\right)\right)$, where ${\mathcal{B}}_j^{\left[1\right]}={\mathcal{H}}_j$ and ${\mathcal{B}}_j^{\left[k\right]}=I-\left(I-{\mathcal{H}}_j\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k-1}}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$ for $k>1$.

3. Update. ${\gamma }_{{\widehat{s}}_k}^{\left[k\right]}={\gamma }_{{\widehat{s}}_k}^{\left[k-1\right]}$ for $j\ne {\widehat{s}}_k$ and ${\gamma }_{{\widehat{s}}_k}^{\left[k\right]}={\gamma }_{{\widehat{s}}_k}^{\left[k-1\right]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_k}^{\left[k\right]}$, where $\nu$ is the step size.

4. Iteration. Repeat step (b)-(c) for $K$ iterations.

5. Stopping. Estimate $\widehat{K}={\text{argmin}}_{1\le k\le K}\text{gMDL}\left({\mathit{RSS}}_{{\widehat{s}}_k}^{\left[k\right]},\text{trace}\left({\mathcal{B}}^{\left[k\right]}\right)\right)$, where ${\mathcal{B}}^{\left[k\right]}=I-\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_k}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$.Thus, ${\gamma }^{\left[\widehat{K}\right]}={\left({\left({\gamma }_0^{\left[\widehat{K}\right]}\right)}^T,\cdots ,{\left({\gamma }_{p-1}^{\left[\widehat{K}\right]}\right)}^T\right)}^T$ is the estimate for $\gamma$ and ${\widehat{\beta }}_j\left(u\right)=B^T\left(u\right){\gamma }_j^{\left[\widehat{K}\right]},j=0,\cdots ,p-1$ are the estimators for varying coefficients. The final estimator of $\tilde{Y}$ is ${\tilde{Y}}^{\left[\widehat{K}\right]}=\tilde{\mathbf{X}}{\gamma }^{\left[\widehat{K}\right]}$.

According to [10] and references therein, the selection of step size $\nu$ is of minor importance as long as it is small. A smaller value of $\nu$ achieves higher prediction accuracy while requires a larger number of boosting iterations and more computing time. A typical value used in literature is $\nu =0.1$.

2.2 Simulation

The performance of the above sparse boosting algorithm is evaluated by studying their performance on simulated data. L2 boosting and sparse boosting methods are compared in their performance of variable selection and function estimation. Sparse boosting method is what we present in this section while L2 boosting method is a relatively simpler version and may not achieve sparse solution in general.

The simulation results from [17] show that both boosting methods can identify important variables while sparse boosting selects much fewer irrelevant variables than L2 boosting. Although in-sample prediction errors (defined as ${\sum }_{i=1}^n{\delta }_i{\left(Y_i-Y_i^{\left[\widehat{K}\right]}\right)}^2/{\sum }_{i=1}^n{\delta }_i$) using L2 boosting is a little bit smaller than using sparse boosting since the former has larger model sizes, the average of root mean integrated squared errors (defined as $\sqrt{\frac{1}{n}{\sum }_{j=0}^5{\sum }_{i=1}^n{\left({\beta }_j\left(u_i\right)-{\widehat{\beta }}_j\left(u_i\right)\right)}^2}$) using sparse boosting is much smaller than that using L2 boosting. Furthermore, when the smoothness assumption in Curry-Schonberg theorem is violated for the coefficient functions, the performance of variable selection remains good. In summary, sparse boosting outperforms L2 boosting in terms of parameter estimation and variable selection.

2.3 Lung cancer data analysis

Lung cancer is the top cancer killer for people in the U.S. Identifying relevant gene expressions in lung cancer is important for treatment and prevention. Our data is from a large multi-site blinded validation study [28] with 442 lung adenocarcinomas. Age is treated as the potential confounder in this analysis, since it is usually strongly correlated with survival time [29]. After removing missing measurements and predictors in overall survival, a total of 439 patients are left in the analysis. For each patient, 22,283 gene expressions are available. The median follow-up time is 46 months (range: 0.03 to 204 months) with the overall censoring rate 46.47 $\%$. The median age at diagnosis is 65 years (range: 33 to 87 years). After adopting a marginal screening procedure to screen out irrelevant genes, variable selection approaches are used to identify important genes associated with lung cancer. With the aim of comparison, except L2 boosting and the proposed sparse boosting, the following existing variable selection approaches for constant-coefficient AFT models are also considered: Buckley-James boosting with linear least squares [23], Buckley-James twin boosting with linear least squares [23], Buckley-James regression with elastic net penalty [30] and SCAD penalty respectively.

The results from [17] show that L2 boosting and sparse boosting for varying-coefficient AFT model not only produce relatively sparser model, but also have smaller in-sample and out-of-sample prediction error compared to the four methods for constant-coefficient AFT model. Again, sparse boosting produce even sparser model than L2 boosting. In conclusion, including age in the varying-coefficient AFT model could lead to more accurate estimate than constant-coefficient AFT model and the proposed sparse boosting method for varying-coefficient AFT model has good performance in terms of estimation, prediction as well as sparsity.

3. Two-step sparse boosting for longitudinal data

Longitudinal data contain repeated measurements collected from the same respondents over time. The assumption that all measurements are independent does not hold for such data. One important question in longitudinal analysis is how to make efficient inference by taking into account of the within subjects correlation. This question has been investigated in depth by many researchers [31, 32] for parametric models. Semiparametric and nonparametric models for longitudinal data are also presented in the literature, see [33, 34]. Recently, there are some development on longitudinal data with high-dimensionalilty using varying-coefficient models [35, 36]. All previous studies adopted the penalty methods.

In this section, a two-step sparse boosting approach is presented to preform the variable selection and the model-based prediction. Specifically, high-dimensional varying-coefficient models with longitudinal data will be considered. In the first step, the sparse boosting approach is utilized to obtain an estimate of the correlation structure. In the second step, the within-subject correlation structure is considered and variable selection and coefficients estimation are achieved by sparse boosting again. The rest of this section is arranged as follows. In Section 3.1, the varying-coefficient model for longitudinal data is formulated and a two-step sparse boosting algorithm is presented. In Section 3.2, simulation studies are conducted to illustrate the validity of the two-step sparse boosting method. In Section 3.3, the performance of two-stage method is assessed by studying yeast cell cycle gene expression data.

3.1 Methodology

3.1.1 Model and estimation

Let $Y_{\mathit{ij}}$ be the continuous outcome for the $j$th measurement of individual $i$ taken at time $t_{\mathit{ij}}\in T$, where $T$ is the time interval on which the measurements are taken. Denote ${\mathbf{X}}_{\mathit{ij}}=\left(X_{\mathit{ij},1},\cdots ,X_{\mathit{ij},p-1}\right)$ to be the corresponding ($p$-$1$)-dimensional covariate vector. The varying-coefficient model which can capture the dynamical impacts of the covariates on the response variable is considered:

$$Y_{\mathit{ij}}={\beta }_0\left(t_{\mathit{ij}}\right)+\sum _{d=1}^{p-1}{X}_{\mathit{ij},d}{\beta }_d\left(t_{\mathit{ij}}\right)+{\epsilon }_{\mathit{ij}},i=1,\cdots ,n,j=1,\cdots ,n_i,$$E10

where ${\beta }_0\left(.\right),{\beta }_1\left(.\right),\cdots ,{\beta }_{p-1}\left(.\right)$ are the unknown smooth coefficient functions of time and ${\epsilon }_i={\left({\epsilon }_{i1},\cdots ,{\epsilon }_{{\mathit{\text{in}}}_i}\right)}^T,i=1,\cdots ,n$ are multivariate error terms with mean zero. Errors are assumed to be uncorrelated for different $i$, but components of ${\epsilon }_i$ are correlated with each other. Without loss of generality, the balanced longitudinal study is considered in the following implementation, i.e., $t_{\mathit{ij}}=t_{\mathit{kj}}$, and $n_i=m$ for all $i$.

The estimation procedure is presented below. In the first step, the within-subject correlation is ignored first and the coefficients are estimated by minimizing the following least squares loss function:

$$\sum _{i=1}^n\sum _{j=1}^m{\left(Y_{\mathit{ij}}-{\beta }_0\left(t_{\mathit{ij}}\right)-\sum _{d=1}^{p-1}{X}_{\mathit{ij},d}{\beta }_d\left(t_{\mathit{ij}}\right)\right)}^2.$$E11

The B-spline basis is used to estimate the coefficient functions ${\beta }_0\left(.\right),{\beta }_1\left(.\right),\cdots ,{\beta }_{p-1}\left(.\right)$. Denote $B\left(.\right)={\left(B_1\left(.\right),\dots ,B_L\left(.\right)\right)}^T$ to be an equal-spaced B-spline basis of dimension $L$. Under certain smoothness assumptions, function ${\beta }_d\left(.\right)$ can be approximated by

$${\beta }_d\left(.\right)\approx B^T\left(.\right){\gamma }_d,d=0,\cdots ,p-1,$$E12

where ${\gamma }_d$ is a loading vector of length $L$. Then the least squares loss function Eq. (11) is close to

$$\sum _{i=1}^n\sum _{j=1}^m{\left(Y_{\mathit{ij}}-B^T\left(t_{\mathit{ij}}\right){\gamma }_0-\sum _{d=1}^{p-1}{X}_{\mathit{ij},d}{B}^T\left(t_{\mathit{ij}}\right){\gamma }_d\right)}^2.$$E13

Further denote $Y_i={\left(Y_{i1},\cdots ,Y_{\mathit{im}}\right)}^T$, $Y={\left(Y_1^T,\cdots ,Y_n^T\right)}^T$, $X_{\mathit{ij},0}=1$, ${\tilde{\mathbf{X}}}_{i,d}={\left(B\left(t_{i1}\right)X_{i1,d},\cdots ,B\left(t_{\mathit{im}}\right)X_{\mathit{im},d}\right)}^T$, ${\tilde{\mathbf{X}}}_i=\left({\tilde{\mathbf{X}}}_{i,0},\cdots ,{\tilde{\mathbf{X}}}_{i,p-1}\right)$, $\tilde{\mathbf{X}}={\left({\tilde{\mathbf{X}}}_1^T,\cdots ,{\tilde{\mathbf{X}}}_n^T\right)}^T$ and $\gamma ={\left({\gamma }_0^T,\cdots ,{\gamma }_{p-1}^T\right)}^T$. Then the target function Eq. (13) can be expressed in the matrix format:

$$\sum _{i=1}^n{\left(Y_i-{\tilde{\mathbf{X}}}_i\gamma \right)}^T\left(Y_i-{\tilde{\mathbf{X}}}_i\gamma \right)\equiv {\left(Y-\tilde{\mathbf{X}}\gamma \right)}^T\left(Y-\tilde{\mathbf{X}}\gamma \right).$$E14

Denote ${\gamma }^{\left[\widehat{K_1}\right]}$ to be the estimator of $\gamma$ by sparse boosting with squared loss function Eq. (14) being loss function, where $\widehat{K_1}$ is the estimated stopping iterations in this step. There is no exact closed form for ${\gamma }^{\left[\widehat{K_1}\right]}$ since it is derived from an iterative algorithm. However it can be evaluated very fast in a computer implementation. The detailed algorithm will be presented in the next subsection.

The first step coefficient estimates are given by

$${\tilde{\beta }}_d\left(t\right)=B^T\left(t\right){\gamma }_d^{\left[\widehat{K_1}\right]},d=0,\cdots ,p-1.$$E15

Write ${\widehat{\epsilon }}_i=Y_i-{\tilde{\mathbf{X}}}_i{\gamma }^{\left[\widehat{K_1}\right]},i=1,\cdots ,n$. The $m\times m$ covariance matrix $\mathit{Cov}\left(Y_i\right)\equiv \mathrm{\Sigma }$ can be estimated by the following empirical estimator

$$\widehat{\mathrm{\sum }}=\frac{1}{n}\sum _{i=1}^n{\widehat{\epsilon }}_i{\widehat{\epsilon }}_i^T.$$E16

In the second step, the estimated correlation structure within repeated measurements is taken into account to form the weighted least squares loss function as follows:

$$\sum _{i=1}^n{\left(Y_i-{\tilde{\mathbf{X}}}_i{\gamma }^{\star }\right)}^T{\widehat{\mathrm{\Sigma }}}^{-1}\left(Y_i-{\tilde{\mathbf{X}}}_i{\gamma }^{\star }\right)\equiv {\left(Y-\tilde{\mathbf{X}}{\gamma }^{\star }\right)}^{T}W\left(Y-\tilde{\mathbf{X}}{\gamma }^{\star }\right),$$E17

where $W=diag\left({\widehat{\mathrm{\Sigma }}}^{-1},\cdots ,{\widehat{\mathrm{\Sigma }}}^{-1}\right)$ is the estimated $\left(n\times m\right)\times \left(n\times m\right)$ weight matrix. Denote ${\gamma }^{\star \left[\widehat{K_2}\right]}$ to be the estimator of ${\gamma }^{\star }$ by sparse boosting with weighted loss function Eq. (17) being the loss function, where $\widehat{K_2}$ is the estimated stopping iterations in the second step. Then the coefficient estimates from the second step are given by

$${\widehat{\beta }}_d\left(t\right)=B^T\left(t\right){\gamma }_d^{\star \left[\widehat{K_2}\right]},d=0,\cdots ,p-1.$$E18

The reliable estimates for the coefficient functions could then be obtained. More details about how to use sparse boosting to get ${\gamma }^{\left[\widehat{K_1}\right]}$ and ${\gamma }^{\star \left[\widehat{K_2}\right]}$ are provided in the following subsection.

3.1.2 Two-step sparse boosting techniques

gMDL can be adopted as the penalized empirical risk function to estimate the update criterion in each iteration and the stopping criterion. gMDL can be expressed in the following form:

$$\begin{array}{ll}\text{gMDL}\left(\mathit{RSS},\text{trace}\left(\mathcal{B}\right)\right)& =log\left(F\right)+\frac{\text{trace}\left(\mathcal{B}\right)}{n\times m}log\left(\frac{Y^{T}Y-\mathit{RSS}}{\text{trace}\left(\mathcal{B}\right)\times F}\right),\\ F& =\frac{\mathit{RSS}}{n\times m-\text{trace}\left(\mathcal{B}\right)},\end{array}$$E19

where $\mathcal{B}$ is the boosting operator and $\mathit{RSS}$ is the residual sum of squares.

The two-step sparse boosting approach is presented more specifically. In the first step, the start value of $\gamma$ is set to zero vector, i.e. ${\gamma }^{\left[0\right]}=\mathbf{0}$, and in each of the $k_1$th iteration ($0<k_1\le K_1$, and $K_1$ is the maximum number of iterations considered in the first step), the residual $R^{\left[k_1\right]}=Y-\tilde{\mathbf{X}}{\gamma }^{\left[k_1-1\right]}$ in present iteration is used to fit each of the $d$th component ${\tilde{\mathbf{X}}}_{,d}={\left({\tilde{\mathbf{X}}}_{1,d}^T,\cdots ,{\tilde{\mathbf{X}}}_{n,d}^T\right)}^T,d=0,\cdots ,p-1$ by treating all the within-subject observations uncorrelated. Then the fit denoted by ${\widehat{\lambda }}_d^{\left[k_1\right]}$ can be calculated by minimizing the squared loss function ${\left(R^{\left[k_1\right]}-{\tilde{\mathbf{X}}}_{,d}\lambda \right)}^T\left(R^{\left[k_1\right]}-{\tilde{\mathbf{X}}}_{,d}\lambda \right)$ with respect to $\lambda$. Therefore, the least squares estimate is ${\widehat{\lambda }}_d^{\left[k_1\right]}={\left[{\left({\tilde{\mathbf{X}}}_{,d}\right)}^T\left({\tilde{\mathbf{X}}}_{,d}\right)\right]}^{-1}{\left({\tilde{\mathbf{X}}}_{,d}\right)}^T{R}^{\left[k_1\right]}$, the corresponding hat matrix is ${\mathcal{H}}_d=\left({\tilde{\mathbf{X}}}_{,d}\right){\left[{\left({\tilde{\mathbf{X}}}_{,d}\right)}^T\left({\tilde{\mathbf{X}}}_{,d}\right)\right]}^{-1}{\left({\tilde{\mathbf{X}}}_{,d}\right)}^T$ and the residual sum of squares is ${\mathit{RSS}}_d^{\left[k_1\right]}={\left(R^{\left[k_1\right]}-{\tilde{\mathbf{X}}}_{,d}{\widehat{\lambda }}_d^{\left[k_1\right]}\right)}^T\left(R^{\left[k_1\right]}-{\tilde{\mathbf{X}}}_{,d}{\widehat{\lambda }}_d^{\left[k_1\right]}\right)$. The chosen element ${\widehat{s}}_{k_1}$ is attained by:

$${\widehat{s}}_{k_1}={\text{argmin}}_{0\le d\le p-1}\text{gMDL}\left({\mathit{RSS}}_d^{\left[k_1\right]},\text{trace}\left({\mathcal{B}}_d^{\left[k_1\right]}\right)\right),$$E20

where ${\mathcal{B}}_d^{\left[1\right]}={\mathcal{H}}_d$ and ${\mathcal{B}}_d^{\left[k_1\right]}=I-\left(I-{\mathcal{H}}_d\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k_1-1}}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$ for $k_1>1$ is the first step boosting operator for choosing $d$th element in the $k_1$th iteration. Hence, there is an unique element ${\tilde{\mathbf{X}}}_{,{\widehat{s}}_{k_1}}$ to be selected at each iteration, and only the corresponding coefficient vector ${\gamma }_{{\widehat{s}}_{k_1}}^{\left[k_1\right]}$ changes, i.e., ${\gamma }_{{\widehat{s}}_{k_1}}^{\left[k_1\right]}={\gamma }_{{\widehat{s}}_{k_1}}^{\left[k_1-1\right]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_{k_1}}^{\left[k_1\right]}$, where $\nu$ is the pre-specified step-size parameter. All the other ${\gamma }_d^{\left[k_1\right]}$ for $d\ne {\widehat{s}}_{k_1}$ keep unchanged. This procedure is repeated for $K_1$ times and the number of iterations $K_1$ can be estimated by

$$\widehat{K_1}={\text{argmin}}_{1\le k_1\le K_1}\text{gMDL}\left({\mathit{RSS}}_{{\widehat{s}}_{k_1}}^{\left[k_1\right]},\text{trace}\left({\mathcal{B}}^{\left[k_1\right]}\right)\right),$$E21

where ${\mathcal{B}}^{\left[k_1\right]}=I-\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k_1}}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$.

From the first step of sparse boosting, the estimator of $\gamma$ is obtained by ${\gamma }^{\left[\widehat{K_1}\right]}={\left({\left({\gamma }_0^{\left[\widehat{K_1}\right]}\right)}^T,\cdots ,{\left({\gamma }_{p-1}^{\left[\widehat{K_1}\right]}\right)}^T\right)}^T$. Then the weight matrix $W$ can be easily obtained too.

In the second step, sparse boosting is used again by taking into account of the correlation structure estimator for the repeated measurements estimated in the first step. The initial value of ${\gamma }^{\star }$ is set to be the coefficient estimator from the first step of sparse boosting, i.e. ${\gamma }^{\star \left[0\right]}={\gamma }^{\left[\widehat{K_1}\right]}$, and in each of the $k_2$th iteration ($0<k_2\le K_2$, and $K_2$ is the maximum number of iterations under consideration in the second step), the residual $R^{\star \left[k_2\right]}=Y-\tilde{\mathbf{X}}{\gamma }^{\star \left[k_2-1\right]}$ in current iteration is used to fit each of the $d$th working element ${\tilde{\mathbf{X}}}_{,d},d=0,\cdots ,p-1$ by incorporating the within-subject correlation estimator from the first step. Then the fit denoted by ${\widehat{\lambda }}_d^{\star \left[k_2\right]}$ can be obtained by minimizing the weighted squared loss function ${\left(R^{\star \left[k_2\right]}-{\tilde{\mathbf{X}}}_{,d}\lambda \right)}^{T}W\left(R^{\star \left[k_2\right]}-{\tilde{\mathbf{X}}}_{,d}\lambda \right)$ with respect to $\lambda$. Thus, the weighted least squares estimate is ${\widehat{\lambda }}_d^{\star \left[k_2\right]}={\left[{\left({\tilde{\mathbf{X}}}_{,d}\right)}^{T}W\left({\tilde{\mathbf{X}}}_{,d}\right)\right]}^{-1}{\left({\tilde{\mathbf{X}}}_{,d}\right)}^T{\mathit{WR}}^{\star \left[k_2\right]}$, the corresponding hat matrix is ${\mathcal{H}}_d^{\star }=\left({\tilde{\mathbf{X}}}_{,d}\right){\left[{\left({\tilde{\mathbf{X}}}_{,d}\right)}^{T}W\left({\tilde{\mathbf{X}}}_{,d}\right)\right]}^{-1}{\left({\tilde{\mathbf{X}}}_{,d}\right)}^{T}W$ and the weighted residual sum of squares is ${\mathit{RSS}}_d^{\star \left[k_2\right]}={\left(R^{\star \left[k_2\right]}-{\tilde{\mathbf{X}}}_{,d}{\widehat{\lambda }}_d^{\star \left[k_2\right]}\right)}^{T}W\left(R^{\star \left[k_2\right]}-{\tilde{\mathbf{X}}}_{,d}{\widehat{\lambda }}_d^{\star \left[k_2\right]}\right)$. The chosen element ${\widehat{s}}_{k_2}$ can be obtained by:

$$\begin{array}{l}{\widehat{s}}_{k_2}={\text{argmin}}_{0\le d\le p-1}\text{gMDL}\left({\mathit{RSS}}_d^{\star \left[k_2\right]},\text{trace}\left({\mathcal{B}}_d^{\star \left[k_2\right]}\right)\right),\end{array}$$E22

where ${\mathcal{B}}_d^{\star \left[1\right]}=I-\left(I-{\mathcal{B}}^{\left[\widehat{K_1}\right]}\right)\left(I-{\mathcal{H}}_d^{\star }\right)$ and ${\mathcal{B}}_d^{\star \left[k_2\right]}=I-\left(I-{\mathcal{B}}^{\left[\widehat{K_1}\right]}\right)\left(I-{\mathcal{H}}_d^{\star }\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k_2-1}}^{\star }\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}^{\star }\right)$ for $k_2>1$ is the second step boosting operator for choosing $d$th element in the $k_2$th iteration. Thus, there is an unique element ${\tilde{\mathbf{X}}}_{,{\widehat{s}}_{k_2}}$ to be selected at each time, and only the corresponding coefficient vector ${\gamma }_{{\widehat{s}}_{k_2}}^{\star \left[k_2\right]}$ change, i.e., ${\gamma }_{{\widehat{s}}_{k_2}}^{\star \left[k_2\right]}={\gamma }_{{\widehat{s}}_{k_2}}^{\star \left[k_2-1\right]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_{k_2}}^{\star \left[k_2\right]}$. While all the other ${\gamma }_d^{\star \left[k_2\right]}$ for $d\ne {\widehat{s}}_{k_2}$ remain the same. This procedure is repeated for $K_2$ times and the estimated stopping iterations $\widehat{K_2}$ is

$$\widehat{K_2}={\text{argmin}}_{1\le k_2\le K_2}\text{gMDL}\left({\mathit{RSS}}_{{\widehat{s}}_{k_2}}^{\star \left[k_2\right]},\text{trace}\left({\mathcal{B}}^{\star \left[k_2\right]}\right)\right),$$E23

where ${\mathcal{B}}^{\star \left[k_2\right]}=I-\left(I-{\mathcal{B}}^{\left[\widehat{K_1}\right]}\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k_2}}^{\star }\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}^{\star }\right)$.

From the second step of sparse boosting, the estimator of ${\gamma }^{\star }$ is arrived by ${\gamma }^{\star \left[\widehat{K_2}\right]}={\left({\left({\gamma }_0^{\star \left[\widehat{K_2}\right]}\right)}^T,\cdots ,{\left({\gamma }_{p-1}^{\star \left[\widehat{K_2}\right]}\right)}^T\right)}^T$. The two-step sparse boosting algorithm for varying-coefficient model with longitudinal data can be summarized in the following form:

Two-step Sparse Boosting Algorithm with Longitudinal Data.

Step I: Use sparse boosting to estimate covariance matrix.

1. Initialization. Let $k_1=0$ and ${\gamma }_0^{\left[k_1\right]}=\mathbf{0},\cdots ,{\gamma }_{p-1}^{\left[k_1\right]}=\mathbf{0}$.

2. Increase $k_1$ by 1. Calculate ${\widehat{s}}_{k_1}={\text{argmin}}_{0\le d\le p-1}\text{gMDL}\left({\mathit{RSS}}_d^{\left[k_1\right]},\text{trace}\left({\mathcal{B}}_d^{\left[k_1\right]}\right)\right)$, where ${\mathcal{B}}_d^{\left[1\right]}={\mathcal{H}}_d$ and ${\mathcal{B}}_d^{\left[k_1\right]}=I-\left(I-{\mathcal{H}}_d\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k_1-1}}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$ for $k_1>1$.

3. Update. ${\gamma }_{{\widehat{s}}_{k_1}}^{\left[k_1\right]}={\gamma }_{{\widehat{s}}_{k_1}}^{\left[k_1-1\right]}$ for $d\ne {\widehat{s}}_{k_1}$ and ${\gamma }_{{\widehat{s}}_{k_1}}^{\left[k_1\right]}={\gamma }_{{\widehat{s}}_{k_1}}^{\left[k_1-1\right]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_{k_1}}^{\left[k_1\right]}$, where $\nu$ is the step-size parameter.

4. Iteration. Repeat step (b)-(c) for some large iteration number $K_1$.

5. Stopping. The optimal iteration number can be taken as $\widehat{K_1}={\text{argmin}}_{1\le k_1\le K_1}\text{gMDL}\left({\mathit{RSS}}_{{\widehat{s}}_{k_1}}^{\left[k_1\right]},\text{trace}\left({\mathcal{B}}^{\left[k_1\right]}\right)\right)$, where ${\mathcal{B}}^{\left[k_1\right]}=I-\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k_1}}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$.

Thus, ${\gamma }^{\left[\widehat{K_1}\right]}={\left({\left({\gamma }_0^{\left[\widehat{K_1}\right]}\right)}^T,\cdots ,{\left({\gamma }_{p-1}^{\left[\widehat{K_1}\right]}\right)}^T\right)}^T$ is the first step estimator for $\gamma$ from sparse boosting and ${\tilde{\beta }}_d\left(t\right)=B^T\left(t\right){\gamma }_d^{\left[\widehat{K_1}\right]}$, $d=0,\cdots ,p-1$ are the varying coefficient estimates ignoring the within-subject correlation. $\mathit{Cov}\left(Y_i\right)$ can be estimated by

$$\widehat{\mathrm{\Sigma }}=\frac{1}{n}\sum _{i=1}^n\left(Y_i-{\tilde{\mathbf{X}}}_i{\gamma }^{\left[\widehat{K_1}\right]}\right){\left(Y_i-{\tilde{\mathbf{X}}}_i{\gamma }^{\left[\widehat{K_1}\right]}\right)}^T.$$

Step II: Use sparse boosting again by incorporating covariance matrix estimator.

1. Initialization. Let $k_2=0$ and ${\gamma }^{\star \left[k_2\right]}={\gamma }^{\left[\widehat{K_1}\right]}$.

2. Increase $k_2$ by 1. Calculate ${\widehat{s}}_{k_2}={\text{argmin}}_{0\le d\le p-1}\text{g MDL}\left({\mathit{RSS}}_d^{\star \left[k_2\right]},\text{trace}\left({\mathcal{B}}_d^{\star \left[k_2\right]}\right)\right)$, where ${\mathcal{B}}_d^{\star \left[1\right]}=I-\left(I-{\mathcal{B}}^{\left[\widehat{K_1}\right]}\right)\left(I-{\mathcal{H}}_d^{\star }\right)$ and ${\mathcal{B}}_d^{\star \left[k_2\right]}=I-\left(I-{\mathcal{B}}^{\left[\widehat{K_1}\right]}\right)\left(I-{\mathcal{H}}_d^{\star }\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k_2-1}}^{\star }\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}^{\star }\right)$ for $k_2>1$.

3. Update. ${\gamma }_{{\widehat{s}}_{k_2}}^{\star \left[k_2\right]}={\gamma }_{{\widehat{s}}_{k_2}}^{\star \left[k_2-1\right]}$ for $d\ne {\widehat{s}}_{k_2}$ and ${\gamma }_{{\widehat{s}}_{k_2}}^{\star \left[k_2\right]}={\gamma }_{{\widehat{s}}_{k_2}}^{\star \left[k_2-1\right]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_{k_2}}^{\star \left[k_2\right]}$.

4. Iteration. Repeat step (b)-(c) for some large iteration number $K_2$.

5. Stopping. The optimal iteration number can be taken as $\widehat{K_2}={\text{argmin}}_{1\le k_2\le K_2}\text{gMDL}\left({\mathit{RSS}}_{{\widehat{s}}_{k_2}}^{\star \left[k_2\right]},\text{trace}\left({\mathcal{B}}^{\star \left[k_2\right]}\right)\right)$, where ${\mathcal{B}}^{\star \left[k_2\right]}=I-\left(I-{\mathcal{B}}^{\left[\widehat{K_1}\right]}\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k_2}}^{\star }\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}^{\star }\right)$.

Therefore, ${\gamma }^{\star \left[\widehat{K_2}\right]}={\left({\left({\gamma }_0^{\star \left[\widehat{K_2}\right]}\right)}^T,\cdots ,{\left({\gamma }_{p-1}^{\star \left[\widehat{K_2}\right]}\right)}^T\right)}^T$ and ${\widehat{\beta }}_d\left(t\right)=B^T\left(t\right){\gamma }_d^{\star \left[\widehat{K_2}\right]}$, $d=0,\cdots ,p-1$ are the final estimator for ${\gamma }^{\star }$ and varying coefficient estimates by the two-step sparse boosting. The final estimate for $Y$ is $\widehat{Y}=\tilde{\mathbf{X}}{\gamma }^{\star \left[\widehat{K_2}\right]}$.

3.2 Simulation

Simulation studies are conducted to evaluate the performance of the above two-step sparse boosting algorithm. The following four methods are compared in terms of variable selection and function estimation performance. M1: two-step L2 boosting (use squared loss for update criterion and gMDL for stopping criterion); M2: two-step sparse boosting; M3: two-step lasso (performs lasso regression in the first step to calculate the estimated within-subject correlation structure using Eq. (14), and use lasso regression in the second step by taking into account of the estimated correlation structure) and M4: two-step elastic net regression (similar as M3 with the elastic net mixing parameter 0.5).

The simulation results from [18] show that all methods are able to identify important variables. However, in terms of sparsity, the two-step sparse boosting method preforms best with smallest number of false positives. Both penalization methods select much more irrelevant variables than boosting methods, with elastic net selects the most. For two-step sparse boosting, results of variable selection are quite stable from step I to step II but for the other approaches, the false positives and thus the sizes of model from step I to step II are expanding. Two-step sparse boosting yields smallest bias for the coefficients estimation among the competing methods. The refined estimates after incorporating the within-subject correlation generally perform better than the initial estimates without taking into account of the within-subject correlation since the two-step methods gain reduction of bias, especially when the within-subject correlation is high. In other words, the reduction of bias from step I to step II are much larger when the within-subject correlation is higher. This is intuitive as in the second step, the within-subject correlation structure estimated from the first step have been taken into account. The similar results obtained for the bias of the estimated covariance matrix. The bias under smaller within-subject correlation is smaller than under larger within-subject correlation. The two-step sparse boosting yields smaller bias of the estimated covariance matrix than other competing methods when the within-subject correlation is high. In summary, the performance of variable selection and functional coefficients estimation for two-step sparse boosting is quite satisfactory.

3.3 Yeast cell cycle gene expression data analysis

The cell cycle is one of the most important activities in life by which cells grow, replicate their chromosomes, undergo mitosis, and split into daughter cells. Thus, identifying cell cycle-regulated genes becomes very important. Adopting a model-based approach, Luan and Li [37] identified $n=297$ cell cycle-regulated genes based on the $\alpha$-factor synchronization experiments. All gene expression levels were measured at $m=18$ different time points covering two cell-cycle periods. Using the same subset of the original data as in [38], a total $p=96$ transcriptional factors (TFs) are included as predictors in the downstream analysis. Wei, Huang and Li [39] proved that the effects of the TFs on gene expression levels are time-dependent. After the independence screening by $l^2$-norm [40] to screen out the irrelevant predictors at first step, several methods can be used to identify the key TFs involved in gene regulation. Except two-step L2 boosting and two-step sparse boosting which take into account of the within-subject correlation in the second step, one-step L2 boosting and one-step sparse boosting which ignore the within-subject correlation are also considered for better comparison. Besides, some two-step penalized approaches are also considered: two-step lasso, two-step adaptive lasso and two-step elastic net (the elastic net mixing parameter 0.5).

The results from [18] show that boosting approaches yield sparser model than the penalized methods. Sparse boosting yields even sparser model and smaller errors in terms of estimation and prediction than L2 boosting. Two-step boosting achieves better performance than one-step boosting with smaller estimation and prediction errors. Two-step sparse boosting method yields the most sparse model, with the smallest in-sample and out-of-sample prediction errors compared to other methods. In terms of the selected TFs, there is a significant overlap between two-step sparse boosting and each of the other methods. In conclusion, the two-step sparse boosting approach performs quite well in terms of variable selection, coefficients estimation and prediction and can provide useful information in identifying the important TFs that take part in the network of regulations.

4. Multi-step sparse boosting for subgroup identification

As personalized medicine is gaining popularity, identification of subgroups of the patients that can gain a higher efficacy from the treatment becomes greatly important. Recently, significant statistical approaches have been proposed to identify subgroups of patients who may be suitable for different treatments. Traditionally, subgroup identification is achieved by parametric partitioning approaches such as Bayesian approaches [41] or classification and regression tree (CART) [42]. Recently, recursive partitioning methods gain popularity since they achieve greater generalizability and efficiency. Such methods include MOB [43], PRIM [44], sequential-BATTing [45] and other non-parametric methods. For a detailed literature review of subgroup identification refer to Lipkovich et al. [46]. In this section, a sparse boosting based subgroup identification method is presented in the context of dense longitudinal data.

In particular, a formal subgroup identification method for high-dimensional dense longitudinal data is presented. It incorporates multi-step sparse boosting into the homogeneous pursuit via change point detection. Firstly, sparse boosting algorithm for individual modeling is first performed to obtain initial estimates. Then, change point detection via binary segmentation is used to identify the subgroup structure of patients. Lastly, the model on each identified subgroups is refitted and again sparse boosting is utilized to remove irrelevant predictors and yield reliable final estimates. The rest of the section is organized as follows. In Section 4.1, the subgroup model is formulated and a detailed method for subgroup identification and estimation is presented. In Section 4.2, the subgroup identification technique is evaluated through simulation studies. In Section 4.3, the feasibility and applicability of the approach is validated by studying a wallaby growth dataset.

4.1 Methodology

4.1.1 Patients model

Denote $Y_{\mathit{it}}$ be the continuous measurement of the $t$th follow-up for patient $i$, where $i=1,\cdots ,n$, $t=1,\cdots ,T_i$. Let ${\mathbf{X}}_{\mathit{it}}=\left(X_{\mathit{it},1},\cdots ,X_{\mathit{it},p}\right)$ be the corresponding $p$-dimensional predictors. Assume $n$ patients are independent. The following longitudinal model for the patients is considered:

$$Y_{\mathit{it}}={\tilde{\beta }}_{i,0}+\sum _{j=1}^p{X}_{\mathit{it},j}{\tilde{\beta }}_{i,j}+{\epsilon }_{\mathit{it}},i=1,\cdots ,n,t=1,\cdots ,T_i.$$E24

where ${\epsilon }_i={\left({\epsilon }_{i1},\cdots ,{\epsilon }_{{\mathit{iT}}_i}\right)}^T,i=1,\cdots ,n$ are multivariate error terms with mean zero. Errors are assumed to be uncorrelated for different $i$, but components of ${\epsilon }_i$ are correlated with each other.

Moreover, the model is further assumed to have the following subgroup structure:

$${\tilde{\beta }}_{i,j}=\left\{\begin{array}{cc}{\beta }_{1,j}& \mathit{\text{when}}i\in {\mathrm{\Omega }}_{1,j}\\ {\beta }_{2,j}& \mathit{\text{when}}i\in {\mathrm{\Omega }}_{2,j}\\ ⋮& ⋮\\ {\beta }_{{\mathcal{N}}_j+1,j}& \mathit{\text{when}}i\in {\mathrm{\Omega }}_{{\mathcal{N}}_j+1,j}\end{array}\right.$$E25

The partition for regression coefficient $\left\{{\tilde{\beta }}_{i,j}:1\le i\le n\right\}$ is $\left\{{\mathrm{\Omega }}_{k,j}:1\le k\le {\mathcal{N}}_j+1\right\}$, which is unknown, and thus there are ${\mathcal{N}}_j+1$ subgroups for the $j$th predictor. All patients are divided into at least ${max}_j\left({\mathcal{N}}_j+1\right)$ and at most ${\prod }_{j=0}^p\left({\mathcal{N}}_j+1\right)$ subgroups by the model. The patients in the same subgroup share a similar relationship between the response and the predictors and have the same set of regression coefficients while different subgroups have different overall relationship between response and covariates. The main aim is to investigate the effects of the predictors on the response for different subgroups.

However, if the number of predictors under consideration is much larger than the number of patients and the number of follow-ups, a serious challenge may arise to estimate regression coefficients. Therefore, instead of adopting traditional methods (eg, MLE), sparse boosting method can be used to estimate the regression coefficients. With this, the dimensionality of features can be reduced and the coefficients of parameters can be obtained simultaneously.

4.1.2 Subgroup identification and estimation

Denote ${\tilde{\beta }}_i={\left({\tilde{\beta }}_{i,0},\cdots ,{\tilde{\beta }}_{i,p}\right)}^T$ and $\tilde{\beta }={\left({\tilde{\beta }}_1^T,\cdots ,{\tilde{\beta }}_n^T\right)}^T$. Firstly, an initial estimator for ${\tilde{\beta }}_i$ is calculated for each subject $i$ through sparse boosting approach using his or her own repeated measurements data; then, homogeneity pursuit via change point detection can be used to identify the change points among ${\beta }_{k,j}$s; lastly, the ${\tilde{\beta }}_i$s can be replaced by the identified subgroup structure, and the final estimator of regression coefficients can be obtained by the sparse boosting algorithm again. The steps for estimating ${\tilde{\beta }}_i$ is outlined as below.

In the first step, individualized modeling via sparse boosting is performed. For each of the $i$th individual, the initial coefficients ${\tilde{\beta }}_i$ can be estimated by minimizing the following least squares loss function:

$$\sum _{t=1}^{T_i}{\left(Y_{\mathit{it}}-{\tilde{\beta }}_{i,0}-\sum _{j=1}^p{X}_{\mathit{it},j}{\tilde{\beta }}_{i,j}\right)}^2.$$E26

Let $Y_i={\left(Y_{i1},\cdots ,Y_{{\mathit{iT}}_i}\right)}^T$, $X_{\mathit{it},0}=1$, $X_{i,j}={\left(X_{i1,j},\cdots ,X_{{\mathit{iT}}_i,j}\right)}^T$, ${\mathbf{X}}_i=\left(X_{i,0},\cdots ,X_{i,p}\right)$. Then the function Eq. (26) can be written in the matrix form:

$${\left(Y_i-{\mathbf{X}}_i{\tilde{\beta }}_i\right)}^T\left(Y_i-{\mathbf{X}}_i{\tilde{\beta }}_i\right).$$E27

Denote ${\tilde{\beta }}_i^{\left[{\widehat{L}}_i\right]}={\left({\tilde{\beta }}_{i,0}^{\left[{\widehat{L}}_i\right]}\cdots ,{\tilde{\beta }}_{i,p}^{\left[{\widehat{L}}_i\right]}\right)}^T$ to be the estimator of ${\tilde{\beta }}_i$ by sparse boosting with Eq. (27) being loss function, where ${\widehat{L}}_i$ is the estimated stopping iterations in this step. This is the initial estimator of ${\tilde{\beta }}_i$. The detailed sparse boosting algorithm will be presented in the next subsection.

In the second step, homogeneity pursuit via change point detection is performed. Binary segmentation algorithm [47] is used to detect the change points among ${\tilde{\beta }}_{i,j}$, $i=1,\cdots ,n$ and to identify the subgroup structure. Let ${\tilde{\beta }}_{i,j}^{\left[{\widehat{L}}_i\right]}$ be the $\left(j+1\right)$th component of ${\tilde{\beta }}_i^{\left[{\widehat{L}}_i\right]}$. For the $j$th covariate, ${\tilde{\beta }}_{i,j}^{\left[{\widehat{L}}_i\right]}$, $i=1,\cdots ,n$, are sorted in ascending order, and denoted by $b_{\left(1\right)}\le \cdots \le b_{\left(n\right)}$. Denote $r_{i,j}$ be the rank of ${\tilde{\beta }}_{i,j}^{\left[{\widehat{L}}_i\right]}$.

For any $1\le l_1<l_2\le n$, denote the scaled difference between the partial means of the first $\tau -l_1+1$ observations and the last $l_2-\tau$ observations to be

$$H_{l_1{l}_2}\left(\tau \right)=\sqrt{\frac{\left(l_2-\tau \right)\left(\tau -l_1+1\right)}{l_2-l_1+1}}\left(\frac{{\sum }_{i=\tau +1}^{l_2}{b}_{\left(l\right)}}{l_2-\tau }-\frac{{\sum }_{i=l_1}^{\tau }{b}_{\left(i\right)}}{\tau -l_1+1}\right).$$E28

Denote $\delta$ to be the threshold, which is a tuning parameter and can be selected by AIC or BIC, then the binary segmentation algorithm is as follows:

1. Find ${\widehat{t}}_1$ such that

   $$H_{1,n}\left({\widehat{t}}_1\right)=\underset{1\le \tau <n}{max}{H}_{1,n}\left(\tau \right).$$E29

   If $H_{1,n}\left({\widehat{t}}_1\right)\le \delta$, there is no change points among $b_{\left(l\right)}$, $l=1,\cdots ,n$, and the change point detection process terminates. Otherwise, ${\widehat{t}}_1$ is added to the set of change points and the region $\left\{\tau :1\le \tau \le n\right\}$ is divided into two subregions: $\left\{\tau :1\le \tau \le {\widehat{t}}_1\right\}$ and $\left\{\tau :{\widehat{t}}_1+1\le \tau \le n\right\}$.

2. Find the change points in the two subregions derived in part (1), respectively. Consider the region $\left\{\tau :1\le \tau \le {\widehat{t}}_1\right\}$ first. Find ${\widehat{t}}_2$ such that

   $$H_{1,{\widehat{t}}_1}\left({\widehat{t}}_2\right)=\underset{1\le \tau <{\widehat{t}}_1}{max}{H}_{1,{\widehat{t}}_1}\left(\tau \right).$$E30

   If $H_{1,{\widehat{t}}_1}\left({\widehat{t}}_2\right)\le \delta$, there is no change point in the region $\left\{\tau :1\le \tau \le {\widehat{t}}_1\right\}$. Otherwise, add ${\widehat{t}}_2$ to the set of change points and divide the region $\left\{\tau :1\le \tau \le {\widehat{t}}_1\right\}$ into two subregions: $\left\{\tau :1\le \tau \le {\widehat{t}}_2\right\}$ and $\left\{\tau :{\widehat{t}}_2+1\le \tau \le {\widehat{t}}_1\right\}$. Similarly, for the region $\left\{\tau :{\widehat{t}}_1+1\le \tau \le n\right\}$, ${\widehat{t}}_3$ can be found such that

   $$H_{{\widehat{t}}_1+1,n}\left({\widehat{t}}_3\right)=\underset{{\widehat{t}}_1+1\le \tau <n}{max}{H}_{{\widehat{t}}_1+1,n}\left(\tau \right).$$E31

   If $H_{{\widehat{t}}_1+1,n}\left({\widehat{t}}_3\right)\le \delta$, there is no change point in the region $\left\{\tau :{\widehat{t}}_1+1\le \tau \le n\right\}$. Otherwise, add ${\widehat{t}}_3$ to the set of change points and divide the region $\left\{\tau :{\widehat{t}}_1+1\le \tau \le n\right\}$ into two subregions: $\left\{\tau :{\widehat{t}}_1+1\le \tau \le {\widehat{t}}_3\right\}$ and $\left\{\tau :{\widehat{t}}_3+1\le \tau \le n\right\}$.

3. For each subregion derived in part (2), the above algorithm is repeated for the subregion $\left\{\tau :1\le \tau \le {\widehat{t}}_1\right\}$ or $\left\{\tau :{\widehat{t}}_1+1\le \tau \le n\right\}$ in part (2) until no change point is detected in any subregions.

The estimated locations for change points are sorted in increasing order and denoted by

$${\widehat{t}}_{\left(1\right)}<{\widehat{t}}_{\left(2\right)}<\cdots <{\widehat{t}}_{\left({\widehat{\mathcal{N}}}_j\right)},$$E32

where ${\widehat{\mathcal{N}}}_j$ is the number of detected change points and could be used to estimate ${\mathcal{N}}_j$. Further denote ${\widehat{t}}_{\left(0\right)}=0$, and ${\widehat{t}}_{\left({\widehat{\mathcal{N}}}_j+1\right)}=n$. Let. ${\widehat{R}}_{i,j}=\left\{\mathrm{\ell }:{\widehat{t}}_{\left(\mathrm{\ell }-1\right)}<r_{i,j}\le {\widehat{t}}_{\left(\mathrm{\ell }\right)}\right\},1\le \mathrm{\ell }\le {\widehat{\mathcal{N}}}_{j+1},$ where $\left\{{\widehat{R}}_{i,j}:1\le i\le n\right\}$ can be used to estimate the grouping index $\left\{R_{i,j}:1\le i\le n\right\}$. The above algorithm can be used to identify the change points for all $j=0,\cdots ,p$ and correspondingly obtain $\left\{{\widehat{R}}_{i,j}:1\le i\le n,0\le j\le p\right\}$. Let $\left\{{\widehat{R}}_{\mathrm{\ell },j}^{\star }:1\le \mathrm{\ell }\le \widehat{\mathcal{N}},0\le j\le p\right\}=\text{unique rows of}\left\{{\widehat{R}}_{i,j}:1\le i\le n,0\le j\le p\right\}$, then $\widehat{\mathcal{N}}$ is the estimated total number of subgroups for patients and the patients index in group $\mathrm{\ell }$ is.

$${\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}=\left\{i:{\widehat{R}}_{i,j}={\widehat{R}}_{\mathrm{\ell },j}^{\star }\right\},1\le \mathrm{\ell }\le \widehat{\mathcal{N}}.$$E33

All the coefficients ${\tilde{\beta }}_{i,j}$s in the same estimated subgroup ${\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}$ are treated to be equal.

In the third step, subgroup modeling is performed by sparse boosting. Incorporating the patients structure identified in step 2, the model is refitted to each of the subgroups via sparse boosting with the following least squares loss function

$$\sum _{i\in {\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}}\sum _{t=1}^{T_i}{\left(Y_{\mathit{it}}-{\tilde{\beta }}_{i,0}-\sum _{j=1}^p{X}_{\mathit{it},j}{\tilde{\beta }}_{i,j}\right)}^2,1\le \mathrm{\ell }\le \widehat{\mathcal{N}}.$$E34

Further denote $Y_{\mathrm{\ell }}^{\star }={\left(Y_{{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\left[1\right]}^T,\cdots ,Y_{{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\left[|{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}|\right]}^T\right)}^T$, $X_{\mathrm{\ell },j}^{\star }={\left(X_{{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\left[1\right],j}^T,\cdots ,X_{{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\left[|{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}|\right],j}^T\right)}^T$, ${\mathbf{X}}_{\mathrm{\ell }}^{\star }=\left(X_{\mathrm{\ell },0}^{\star },\cdots ,X_{\mathrm{\ell },p}^{\star }\right)$ and ${\tilde{\beta }}_{\mathrm{\ell }}^{\star }={\left({\tilde{\beta }}_{{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\left[1\right]}^T,\cdots ,{\tilde{\beta }}_{{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\left[|{\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}|\right]}^T\right)}^T$ for $\mathrm{\ell }=1,\cdots ,\widehat{\mathcal{N}}$, where ${\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\left[i\right]$ is the $i$th element of ${\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}$ and $\mid {\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\mid$ is the number of elements in ${\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}$. The function Eq. (34) can be written in the matrix form:

$${\left(Y_{\mathrm{\ell }}^{\star }-{\mathbf{X}}_{\mathrm{\ell }}^{\star }{\tilde{\beta }}_{\mathrm{\ell }}^{\star }\right)}^T\left(Y_{\mathrm{\ell }}^{\star }-{\mathbf{X}}_{\mathrm{\ell }}^{\star }{\tilde{\beta }}_{\mathrm{\ell }}^{\star }\right),1\le \mathrm{\ell }\le \widehat{\mathcal{N}}.$$E35

Denote ${\tilde{\beta }}_{\mathrm{\ell }}^{\star \left[{\widehat{L}}_{\mathrm{\ell }}^{\star }\right]}$ to be the estimate for ${\tilde{\beta }}_{\mathrm{\ell }}^{\star }$ by sparse boosting with Eq. (35) being the loss function, where ${\widehat{L}}_{\mathrm{\ell }}^{\star }$ is the estimated number of stopping iterations in this step. The estimator for coefficient ${\tilde{\beta }}_i$ is

$${\widehat{\beta }}_i=\left\{{\tilde{\beta }}_{\mathrm{\ell }}^{\star \left[{\widehat{L}}_{\mathrm{\ell }}^{\star }\right]}\text{for}i\in {\widehat{\mathrm{\Omega }}}_{\mathrm{\ell }}\right\},1\le i\le n.$$E36

More details about how to use sparse boosting to obtain $\left\{{\tilde{\beta }}_i^{\left[{\widehat{L}}_i\right]},1\le i\le n\right\}$ and $\left\{{\tilde{\beta }}_{\mathrm{\ell }}^{\star \left[{\widehat{L}}_{\mathrm{\ell }}^{\star }\right]},1\le \mathrm{\ell }\le \widehat{\mathcal{N}}\right\}$ are given in the following subsection.

4.1.3 Multi-step sparse boosting techniques

gMDL can be used as the penalized empirical risk function to estimate the update criterion in each iteration and the stopping criterion to avoid the selection of the tuning parameter. gMDL can be expressed in the following form:

$$\begin{array}{ll}\text{gMDL}\left(Y,\mathit{RSS},\text{trace}\left(\mathcal{B}\right)\right)& =log\left(F\right)+\frac{\text{trace}\left(\mathcal{B}\right)}{\mid Y\mid }log\left(\frac{Y^{T}Y-\mathit{RSS}}{\text{trace}\left(\mathcal{B}\right)\times F}\right),\\ F& =\frac{\mathit{RSS}}{\mid Y\mid -\text{trace}\left(\mathcal{B}\right)},\end{array}$$E37

where $Y$ is the vector of response variable, $\mid Y\mid$ is the length of $Y$, $\mathcal{B}$ is the boosting operator and $\mathit{RSS}$ is the residual sum of squares.

The sparse boosting procedure is described in details. The starting value of ${\tilde{\beta }}_i$ is set to zero vector, i.e. ${\tilde{\beta }}_i^{\left[0\right]}=0$, and in each of the $l_i$th iteration ($0<l_i\le L_i$, and $L_i$ is the maximum number of iterations considered in this step), the residual $R^{\left[l_i\right]}=Y_i-{\mathbf{X}}_i{\tilde{\beta }}_i^{\left[l_i-1\right]}$ in present iteration is used to fit each of the $j$th element $X_{i,j},j=0,\cdots ,p$. The fit denoted by ${\widehat{\lambda }}_j^{\left[l_i\right]}$ can be obtained by minimizing the squared loss function ${\left(R^{\left[l_i\right]}-X_{i,j}\lambda \right)}^T\left(R^{\left[l_i\right]}-X_{i,j}\lambda \right)$ with respect to $\lambda$. Thus, the least squares estimate is ${\widehat{\lambda }}_j^{\left[l_i\right]}={\left[{\left(X_{i,j}\right)}^T\left(X_{i,j}\right)\right]}^{-1}{\left(X_{i,j}\right)}^T{R}^{\left[l_i\right]}$, the corresponding hat matrix is ${\mathcal{H}}_j=\left(X_{i,j}\right){\left[{\left(X_{i,j}\right)}^T\left(X_{i,j}\right)\right]}^{-1}{\left(X_{i,j}\right)}^T$ and the residual sum of squares is ${\mathit{RSS}}_j^{\left[l_i\right]}={\left(R^{\left[l_i\right]}-X_{i,j}{\widehat{\lambda }}_j^{\left[l_i\right]}\right)}^T\left(R^{\left[l_i\right]}-X_{i,j}{\widehat{\lambda }}_j^{\left[l_i\right]}\right)$. The selected entry ${\widehat{s}}_{l_i}$ is obtained by:

$$\begin{array}{l}{\widehat{s}}_{l_i}={\text{argmin}}_{0\le j\le p}\text{gMDL}\left(Y_i,{\mathit{RSS}}_j^{\left[l_i\right]},\text{trace}\left({\mathcal{B}}_j^{\left[l_i\right]}\right)\right),\end{array}$$E38

where ${\mathcal{B}}_j^{\left[1\right]}={\mathcal{H}}_j$ and ${\mathcal{B}}_j^{\left[l_i\right]}=I-\left(I-{\mathcal{H}}_j\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{l_i-1}}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$ for $l_i>1$ is the boosting operator for choosing $j$th entry in the $l_i$th iteration in this step. Hence, there is an unique element $X_{i,{\widehat{s}}_{l_i}}$ to be selected at each iteration, and only the corresponding coefficient vector ${\tilde{\beta }}_{i,{\widehat{s}}_{l_i}}^{\left[l_i\right]}$ changes, i.e., ${\tilde{\beta }}_{i,{\widehat{s}}_{l_i}}^{\left[l_i\right]}={\tilde{\beta }}_{i,{\widehat{s}}_{l_i}}^{\left[l_i-1\right]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_{l_i}}^{\left[l_i\right]}$, where $\nu$ is the pre-specified step-size parameter. All the other ${\tilde{\beta }}_{i,j}^{\left[l_i\right]}$ for $j\ne {\widehat{s}}_{l_i}$ keep unchanged. This procedure is repeated for $L_i$ times and the number of iterations $L_i$ can be estimated by

$${\widehat{L}}_i={\text{argmin}}_{1\le l_i\le L_i}\text{gMDL}\left(Y_i,{\mathit{RSS}}_{{\widehat{s}}_{l_i}}^{\left[l_i\right]},\text{trace}\left({\mathcal{B}}^{\left[l_i\right]}\right)\right),$$E39

where ${\mathcal{B}}^{\left[l_i\right]}=I-\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{l_i}}\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}\right)$.

From the above sparse boosting approach, the estimator of ${\tilde{\beta }}_i$ is ${\tilde{\beta }}_i^{\left[{\widehat{L}}_i\right]}={\left({\tilde{\beta }}_{i,0}^{\left[{\widehat{L}}_i\right]},\cdots ,{\tilde{\beta }}_{i,p}^{\left[{\widehat{L}}_i\right]}\right)}^T$, $i=1,\cdots ,n$. Then the subgroup structure can be obtained by homogeneity pursuit via change point detection.

Next, sparse boosting is used again for each estimated subgroups. The starting value of ${\tilde{\beta }}_{\mathrm{\ell }}^{\star }$ is set to zero vector, i.e. ${\tilde{\beta }}_{\mathrm{\ell }}^{\star \left[0\right]}=0$, and in each of the $l_{\mathrm{\ell }}^{\star }$th iteration ($0<l_{\mathrm{\ell }}^{\star }\le L_{\mathrm{\ell }}^{\star }$, and $L_{\mathrm{\ell }}^{\star }$ is the maximum number of iterations considered in this stage), the residual $R^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}=Y_{\mathrm{\ell }}^{\star }-{\mathbf{X}}_{\mathrm{\ell }}^{\star }{\tilde{\beta }}_{\mathrm{\ell }}^{\star \left[l_{\mathrm{\ell }}^{\star }-1\right]}$ in present iteration is used to fit each of the $j$th component $X_{\mathrm{\ell },j}^{\star },j=0,\cdots ,p$. Then the fit denoted by ${\widehat{\lambda }}_j^{\star \left[l_i^{\star }\right]}$ can be calculated by minimizing the squared loss function ${\left(R^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}-X_{\mathrm{\ell },j}^{\star }\lambda \right)}^T\left(R^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}-X_{i,j}^{\star }\lambda \right)$ with respect to $\lambda$. Therefore, the least squares estimate is ${\widehat{\lambda }}_j^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}={\left[{\left(X_{\mathrm{\ell },j}^{\star }\right)}^T\left(X_{\mathrm{\ell },j}^{\star }\right)\right]}^{-1}{\left(X_{\mathrm{\ell },j}^{\star }\right)}^T{R}^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}$, the corresponding hat matrix is ${\mathcal{H}}_j^{\star }=\left(X_{\mathrm{\ell },j}^{\star }\right){\left[{\left(X_{\mathrm{\ell },j}^{\star }\right)}^T\left(X_{\mathrm{\ell },j}^{\star }\right)\right]}^{-1}{\left(X_{\mathrm{\ell },j}^{\star }\right)}^T$ and the residual sum of squares is ${\mathit{RSS}}_j^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}={\left(R^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}-X_{\mathrm{\ell },j}^{\star }{\widehat{\lambda }}_j^{\left[l_{\mathrm{\ell }}^{\star }\right]}\right)}^T\left(R^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}-X_{\mathrm{\ell },j}^{\star }{\widehat{\lambda }}_j^{\left[l_{\mathrm{\ell }}^{\star }\right]}\right)$. The chosen element ${\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}^{\star }$ is attained by:

$$\begin{array}{l}{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}^{\star }={\text{argmin}}_{0\le j\le p}\text{gMDL}\left(Y_{\mathrm{\ell }}^{\star },{\mathit{RSS}}_j^{\star \left[l_{\mathrm{\ell }}^{\star }\right]},\text{trace}\left({\mathcal{B}}_j^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}\right)\right),\end{array}$$E40

where ${\mathcal{B}}_j^{\star \left[1\right]}={\mathcal{H}}_j^{\star }$ and ${\mathcal{B}}_j^{\star \left[l_{\mathrm{\ell }}\right]}=I-\left(I-{\mathcal{H}}_j^{\star }\right)\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }-1}}^{\star }\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}^{\star }\right)$ for $l_{\mathrm{\ell }}^{\star }>1$ is the boosting operator for choosing $j$th element in the $l_{\mathrm{\ell }}^{\star }$th iteration in this stage. Hence, there is an unique element $X_{\mathrm{\ell },{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}}^{\star }$ to be selected at each iteration, and only the corresponding coefficient vector ${\tilde{\beta }}_{\mathrm{\ell },{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}}^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}$ changes, i.e., ${\tilde{\beta }}_{\mathrm{\ell },{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}}^{\left[l_{\mathrm{\ell }}\right]}={\tilde{\beta }}_{\mathrm{\ell },{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}}^{\star \left[l_{\mathrm{\ell }}^{\star }-1\right]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}}^{\left[l_{\mathrm{\ell }}^{\star }\right]}$, where $\nu$ is the pre-specified step-size parameter. All the other ${\tilde{\beta }}_{\mathrm{\ell },j}^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}$ for $j\ne {\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}$ keep unchanged. This procedure is repeated for $L_{\mathrm{\ell }}^{\star }$ times and the number of iterations $L_{\mathrm{\ell }}^{\star }$ can be estimated by

$${\widehat{L}}_i^{\star }={\text{argmin}}_{1\le l_{\mathrm{\ell }}\le L_i^{\star }}\text{gMDL}\left(Y_{\mathrm{\ell }}^{\star },{\mathit{RSS}}_{{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}}^{\star \left[l_{\mathrm{\ell }}^{\star }\right]},\text{trace}\left({\mathcal{B}}^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}\right)\right),$$E41

where ${\mathcal{B}}^{\star \left[l_{\mathrm{\ell }}^{\star }\right]}=I-\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_{l_{\mathrm{\ell }}^{\star }}}^{\star }\right).\cdots .\left(I-\nu {\mathcal{H}}_{{\widehat{s}}_1}^{\star }\right)$.

From the second step of sparse boosting, the estimator of ${\tilde{\beta }}_{\mathrm{\ell }}$ is ${\tilde{\beta }}_{\mathrm{\ell }}^{\star \left[{\widehat{L}}_{\mathrm{\ell }}^{\star }\right]}={\left({\tilde{\beta }}_{\mathrm{\ell },0}^{\star \left[{\widehat{L}}_{\mathrm{\ell }}^{\star }\right]},\cdots ,{\tilde{\beta }}_{\mathrm{\ell },p}^{\star \left[{\widehat{L}}_{\mathrm{\ell }}^{\star }\right]}\right)}^T$, $\mathrm{\ell }=1,\cdots ,\widehat{\mathcal{N}}$.

4.2 Simulation

Extensive simulations are conducted to evaluate the performance of the proposed procedure. The accuracy of subgrouping, feature selection, coefficients estimation and prediction are assessed in the setting of different number of patients and repeated measurements. To understand the advantage of the proposed method better, the following four approaches are also considered. M1: the homogeneous model fitting method which treats all patients as one group and use sparse boosting for the single model to estimate $\tilde{\beta }$; M2: the heterogeneous model fitting method which uses initial pre-grouping estimate ${\tilde{\beta }}_i^{\left[{\widehat{L}}_i\right]}$ as the final estimate of ${\tilde{\beta }}_i$; M3: same as the proposed method but in step 2, instead of detecting the change points for coefficients of each covariate ${\tilde{\beta }}_{i,j}^{\left[{\widehat{L}}_i\right]}$, $i=1,\cdots ,n$ for $j=0,\cdots ,p$, it detects the change points among ${\left({\tilde{\beta }}_1^{T\left[{\widehat{L}}_1\right]},\cdots ,{\tilde{\beta }}_n^{T\left[{\widehat{L}}_n\right]}\right)}^T$ similarly to Ke et al. [48]; M4: the proposed method.

The results from [19] show that the naive homogeneous model fitting method M1 can rarely identify the important covariates while the over-parameterized model fitting method M2 and other two methods (M3 & M4) which identify subgroup structures consistently yield true positives equal to the true number of important covariates. Compared these three methods which can identify the important covariates, the proposed method produces smallest false positives. In addition, the number of false positives is decreasing when there is an increase in cluster size. Neither the homogeneous model fitting method nor heterogeneous model fitting method is able to identify the true structure among patients. The method M3 produces much more subgroups than it really has, while the proposed method M4 identified the number of subgroups closest to the actual number of subgroups. Furthermore, the probability of identifying the true subgroups becomes larger when the number of repeated measurements increases. For in-sample prediction, the over-parameterized model M2 performs the best while the methods M3 & M4 performs very competitively. However, for out-of-sample prediction, method M4 is the best. M1 is inferior to M4, yielding poor results of estimation and prediction. In summary, the proposed method preforms pretty well in terms of subgroup identification, variable selection, estimation as well as prediction.

4.3 Wallaby growth data analysis

The proposed subgroup identification method is applied to wallaby growth data, which is from the Australian Dataset and Story Library (OzDASL) and can be found at http://www.statsci.org/data/oz/wallaby.html. The data set has 77 Tammar wallabies’ growth measurements which were taken longitudinally. The response variable is the weight of wallabies (tenths of a gram). The predictors involve length of head, ear, body, arm, leg, tail, foot and their second order interactions. Therefore, a total of 35 predictors are included in the analysis. After removing the missing data, 43 Tammar wallabies are kept in our dataset. The number of repeated measurements ranges from 9 to 34 (median: 23). To have a better understanding of the wallabies’ growth trend, the questions of which parts of body would affect the weight and whether the length of each body parts have the same effects on the weight for all wallabies are investigated, i.e. is there any subgroups among wallabies. Except the above subgroup identification method (SB-CPD1), the other 3 methods studied in simulation are also considered, i.e. homogeneous model fitting method (SB-Homogeneous), heterogeneous model fitting method (SB-Heterogeneous) and the method similar to SB-CPD1 but identifying subgroups via other method in Ke et al. [48] (SB-CPD2). In addition, the following subgroup identification methods incorporating penalized methods are also investigated: similar to our proposed method but instead of using sparse boosting, lasso (Lasso-CPD1), elastic net (ElasticNet-CPD1), SCAD (SCAD-CPD1) or MCP (MCP-CPD1) is used.

The results from [19] show that although Lasso-CPD1 and ElasticNet-CPD1 yield smaller in-sample prediction error by keeping all 35 covariates, they have relatively large out-of-sample prediction errors due to over-fitting problem. The subgroup identification method via sparse boosting keeps smaller number of predictors, achieves sparser model than penalized methods. The proposed method SB-CPD1 identifies smaller number of subgroups and predictors than alternative competing methods while produces smallest out-of-sample prediction errors. In conclusion, the proposed subgroup identification method provides a more precise definition for various subgroups. It may also result in a more accurate medical decision making for these subjects.

5. Conclusions

In this chapter, we discussed various sparse boosting based machine learning methods in the context of high-dimensional data problems. Specifically, we presented the sparse boosting procedure and two-step sparse boosting procedure for nonparametric varying-coefficient models with survival data and repeatedly measured longitudinal data respectively to simultaneously perform variable selection and estimation of functional coefficients. We further presented the multi-step sparse boosting based subgroup identification method with longitudinal patient data to identify subgroups that exhibit different treatment effects. The extensive numerical studies show the validity and effectiveness of our proposed methods and the real data analysis further demonstrate their usefulness and advantages.