Improved mixed model for longitudinal data analysis using shrinkage method

The problem of multicollinearity among predictor variables is a frequent issue in longitudinal data analysis. In this context, this paper proposes a mixed ridge regression model via shrinkage methods to analyze such data. Furthermore, in view of obtaining more efficient estimators, we propose preliminary and Stein-type estimators using prior information for fixed-effects parameters. The model parameters are estimated via the EM algorithm. A simulation study is also presented to assess the performance of the estimators under different estimation methods. An application to the HIV data is also illustrated.


Introduction
In longitudinal data setup, repeated measures of some variables of interest are collected over a specified time period for different independent subjects or individuals. Such types of data are commonly encountered in medical research where the responses are subject to various time-dependent and time-constant effects such as pre-and post-treatment types, gender effect, baseline measures and among others (see Mamode Khan et al. [1], Yuan et al. [2], Verbeke et al. [3], Temesgen and Kebede [4], Seyoum et al. [5] and the references therein). It is quite natural, in the above examples, the repeated measures shall exhibit some forms of dependence that may be resulted from some serial or random effects as outlined by Zeger and Liang [6], Thall and Vail [7], Laird and Ware [8], Sutradhar [9] and Sutradhar and Jowaheer [10]. Thus, the main purpose of the longitudinal studies is to estimate the effects of the various parameters and determine their significance while the dependence estimate is treated as secondary. In this context, FitzMaurice and Laird [11] and Sutradhar et al. [12] have proposed various likelihood-based and pseudo-likelihood-based estimation procedures to estimate the regression effects but the efficiency of the estimators in these approaches may be questionable, in particular, under multi-collinearity among the predictor variables as considered by Eliot et al. [13], Hossain et al. [14] and Saleh et al. [15].
Since longitudinal data mostly arise from clinical studies, the expert knowledge about the parameters has vital impact on the output and is thus an important component in the estimation of model parameters. The preliminary test and shrinkage techniques are mostly used mechanisms in which a prior knowledge can be included in the estimation stage (see papers by Ali and Saleh [16], Ahmed and Fallahpour [17] ,Roozbeh and Arashi [18] and Yuzbasi and Ahmed [19], Yuzbasi et al. [20] and Asar [21] and the references therein). In this paper, we develop the preliminary test and shrinkage estimation methods for the analysis of longitudinal data in ridge regression context, where some parameters are subject to certain/uncertain restrictions. By this, we improve the 1 3 estimation technique, in both the mean squared error (MSE) and mean prediction error (MPE) senses.
We begin with the linear mixed effects (LM) model given by w h e r e i = 1, … , n r e p r e s e n t s i n d i v i d u a l , i = (y i1 , y i2 , … , y in i ) T is a vector of n i observations for individual i, and i is the corresponding n i × p design matrix of fixed-effect covariates. We further assume i ∼ N q (0, ) is person-specific random effects, i the corresponding random effects design matrix and i ∼ N n (0, 2 n i ×n i ) . Let , and be appropriately defined matrices representing the concatenation of the corresponding variables over all individuals i. Then, the LM model in matrix notation has form where : N × 1 , : N × p , : p × 1 , : N × nq and : nq × 1 . Here, N = ∑ n i=1 n i . The log-likelihood function of based on this model is given by where = Var( ) = T + 2 and i is the variance component corresponding to individual i. Maximizing (3) with respect to the fixed-effects parameter vector, , in the non-penalized setting is equivalent to minimizing the least squares objective function that gives the estimate of as Consider a situation in which the multi-collinearity problem is present. The ridge regression approach designed specifically to handle correlated predictors involves introducing a shrinkage penalty k to the least squares equation, and subsequently solving for the value of such that The structure of the paper is therefore organized as follows: In "Shrinkage approach" section, the novel estimation approach based on shrinkage method is proposed which also includes some computational aspects and the EM algorithm. "Simulation study" section provides the Monte Carlo simulation experiments followed by the section on the application of the new techniques to the HIV data, and a comparative study is also presented. The conclusion is presented in final section.

Shrinkage approach
In this section, we consider the mixed ridge (MR) model (2) and develop preliminary test and Stein-type estimation of the fixed-effects parameter when it is a priori suspected that = 0 . Usually, comes from different sources: (1) a fact known from theoretical or experimental considerations, (2) hypothesis that may have to be tested or (3) an artificially imposed condition to reduce redundancy in the description of the model. Here, we interpret as expert knowledge.

Incorporating expert knowledge
In the statistical literature, preliminary test estimation of parameters was introduced by Bancroft [22] to estimate the parameters of a model when it is suspected that some "uncertain prior information" (UPI) on the parameter of interest is available. The method involves a statistical test of the UPI based on an appropriate statistic and a decision on whether the model-based sample estimate or the prior information-based estimate of the model parameters should be taken.
In our case, if we suspect to be 0 , then 0 is termed as the restricted estimator (RE) that must be incorporated in the analysis. Otherwise, ̂ MR is used. In practice, the prior information that = 0 is uncertain. The doubt on this prior information can be removed by testing the following hypothesis: As a result of this test, we choose 0 or not based on the rejection or acceptance of 0 . Accordingly, we develop the estimator called the preliminary test MR (PTMR) estimator. Using indicator function, this estimator can be rewritten as where I A is the indicator function of the set A, L is the Wald statistic given by L = (̂ − 0 ) T ( T −1 )(̂ − 0 ) and L n ( ) is the upper level critical value of the -distribution with p d.f. The PTMR estimator is highly dependent on the level of significance and has discrete nature which is simplified to one of the extremes ̂ MR and 0 , according to the output of test. In this respect, making use of a continuous and -free estimator may make more sense.
Stein-type estimation was introduced by Stein [23] and James and Stein [24] in the statistical literature. It combines UPI on the parameters of interest and the sample observation from the statistical model. In the context of LM model, using the same approach as in Saleh [15], the Stein-type estimator of has form Unlike the PTMR estimator, the Stein-type estimator is a smooth function and independent of the level of significance .
Notice that the forms of ̂ S MR and ̂ PT MR are similar where Eq. (5) obtained from Eq. (4) by replacing I(L n < L n ( )) by cL −1 n to make it independent of .

Computational aspects
Consider the setting in which the variance parameters = ( , 2 ) are unknown. Eliot et al. [13] proposed an extension of the expectation-maximization (EM) algorithm described by Laird and Ware [8] that includes an additional step for estimation of the ridge component. We refine their procedure to evaluate shrinkage estimators as well.
Further, for the ease of computations, we use the estimate of Hoerl and Kenard [25] for the ridge parameter as the initial value. It is given by In what follows, we propose the refined EM algorithm.
MR and the sufficient statisticst n i and n is the number of individuals in our sample, and letV (t+1 Steps (i)-(iv) a large number of times and until a convergence criterion is met and obtainβ M andβ MR v) Compute L n usingβ M and

Simulation study
A Monte Carlo simulation study is conducted to evaluate the performance of the proposed PTMR and shrinkage estimators compared to the MR estimator of Eliot et al. [13]. In our simulation scheme, we fix n i = 4 measurements for each subject i and generate data according to the model of (2) where . To generate multicollinear data, using the "EnvStats" package in R, each predictor variable is assumed to arise from N(5, 1), and the correlation between predictor variables is taken from the set {0.0, 0.  Tables 1 and 2. From Table 2, it is clear the estimation is dependent on the level of significance .
From Table 1, it is apparent the shrinkage mixed ridge (shrinkage MR) estimator has smaller MSE and standard error (sd). Hence, the shrinkage MR is the best among all other competitors; i.e., the shrinkage MR performs better than the mixed, MR and shrinkage mixed (shrinkage M). Knowing this, the preliminary test approach is only applied to the mixed ridge estimator, giving rise to the PTMR estimator. According to the results of Table 2, as the level of significance increases, the MSE increases. The graphs of the MSE against the different values of are also shown in Fig. 1.

3
Although as level of multi-collinearity increases, so does the MSE values, the proposed PTMR estimator has smaller MSE among all. Further, the PTMR and shrinkage MR estimators perform better than the M and MR estimators in multi-collinear situations.

HIV data analysis
In a similar framework as explained in Temesgen and Kebede [4], this section focuses on analyzing HIV data using the linear mixed model. In particular, in this study, we analyze the performance of the proposed estimators using the aids dataset taken from "JMbayes" package in R. The dataset consists of seven covariates for each n = 467 patients. The response variable is the CD4, and we consider the gender and prevOI variables as the random effects, in our study. Didanosine versus zalcitabine in HIV patients, a randomized clinical trial data were collected to compare the efficacy and safety of two antiretroviral drugs in treating patients who had failed or were intolerant of zidovudine (AZT) therapy. Table 3 introduces variables of this study: Variables have been measured n i time for each individuals, so we have a longitudinal dataset; some of the variables like gender in this set will put the subjects in special groups, so we can consider these variables as random effects and we can use mixed models to analyze this dataset, but we know the high degree of correlation among predictors is expected, because variables are measured several times, so we have multi-collinearity data; to combat this difficulty, we can use mixed ridge regression (see Eliot et al. [13]). Information on the important variables is summarized in Table 4.
In addition, we use shrinkage methods to increase estimation efficiency. For the purpose of utilizing the mixed model in "Shrinkage approach" section, the log transform was applied to the CD 4 counts.    Table 6, it is clear that the drug ddI is more effective in curing the infected cells than ddC. In fact, the rate of improvement through using ddI increases by 63.41% as compared to using ddC. This conforms with the negative sign of the parameter "obstime" which indicates that the proliferation of the HIV virus is under control as the time point increases by an approximate rate of 5.715% . The negative sign in the "death" parameter also denotes that using ddC, the rate of deaths and the time to survival improve by 47.83% and 8.36% respectively. Besides, the standard deviation estimates illustrate that these explanatory variables are significant, and hence based on the above dataset, the drug ddI (didanosine) is a recommendable remedy.
From the medical point of view as well, it is shown that ddI yields better treatment in controlling the growth of the HIV virus in the human body (see Molina et al. [26]) while the drug ddC (zalcitabine) has been strongly recommended to be unused due to its countereffects as discussed in the book "clinical neurotoxicology" by Dobbs [27] and Bilgrami and O'Keefe [28].   To compare the performance of the shrinkage MR estimator, we evaluate the MPE; the lesser, the better. In what follows, we describe the scheme we used to derive the MPE. For our purpose, a K-fold cross-validation is used to obtain an estimate of the prediction errors of the model. In a K-fold cross-validation, the dataset is randomly divided into K subsets of roughly equal size. One subset is left aside, {( test , test )} , termed as test set, while the remaining K − 1 subsets, called the training set, are used to fit model. The resultant estimator is called ̂ train . The fitted model is then used to predict the responses of test dataset. Finally, prediction errors are obtained by taking the squared deviation of the observed and predicted values in the test set, i.e., The process is repeated for all K subsets, and the prediction errors are combined. To account for the random variation of the cross-validation, the process is reiterated N times and the average prediction error is estimated which is given by where PE k i is the prediction error of considering kth test set in ith iteration. If the value of MPE is lesser, the estimator is preferred, comparatively.
Our results are based on N = 500 case re-sampled bootstrap sample. In Table 5, we report the estimates and MPE  values. Based on the results, the proposed shrinkage MR estimator performs better than the others, in MPE sense. Further, the absolute value of estimates in the shrinkage MR estimates is lesser than the others. The box plots of the PE are shown in Fig. 2. From the results in the estimation table, it could be deduced that the didanosine drug provides a better treatment.

Conclusion
In this paper, we developed a preliminary test and Stein-type ridge regression estimation in linear mixed model for longitudinal data analysis. Hence, we considered a penalized likelihood approach and proposed the shrinkage mixed ridge estimator for the vector of regression coefficients. An EM algorithm is also exhibited to solve the penalized likelihood for the unknown parameters. Simulation studies demonstrated the good performance of the proposed estimator for multicollinear situations compared to the maximum likelihood estimator. In addition, the above model has contributed largely to justify the use of didanosine in improving the health states of HIV patients, as stated in various biomedical studies. Henceforth, such model and its estimation step based on shrinkage is highly commendable for medical studies of such genre.