Improved mixed model for longitudinal data analysis using shrinkage method

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. [ ¹ ], Yuan et al. [ ² ], Verbeke et al. [ ³ ], Temesgen and Kebede [ ⁴ ], Seyoum et al. [ ⁵ ] 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 [ ⁶ ], Thall and Vail [ ⁷ ], Laird and Ware [ ⁸ ], Sutradhar [ ⁹ ] and Sutradhar and Jowaheer [ ¹⁰ ]. 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 [ ¹¹ ] and Sutradhar et al. [ ¹² ] 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. [ ¹³ ], Hossain et al. [ ¹⁴ ] and Saleh et al. [ ¹⁵ ].

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 [ ¹⁶ ], Ahmed and Fallahpour [ ¹⁷ ] ,Roozbeh and Arashi [ ¹⁸ ] and Yuzbasi and Ahmed [ ¹⁹ ], Yuzbasi et al. [ ²⁰ ] and Asar [ ²¹ ] 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 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

The log-likelihood function of Y based on this model is given by

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

Shrinkage approach

In this section, we consider the mixed ridge (MR) model ( ² ) 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.

Computational aspects

Consider the setting in which the variance parameters θ=(D,σ2) are unknown. Eliot et al. [ ¹³ ] proposed an extension of the expectation–maximization (EM) algorithm described by Laird and Ware [ ⁸ ] 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 [ ²⁵ ] for the ridge parameter as the initial value. It is given by

k^0=σ^2α^max2,

where α^max=max(α^1,…,α^p) in which α^=ΓTβ^ and Γ is the orthogonal eigenvector is spectral decomposition of XTX , i.e., XTX=ΓΛΓT , Λ=diag(λ1,…,λp) , where λi is the i th eigenvalue of XTX . In what follows, we propose the refined EM algorithm.

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. [ ¹³ ]. In our simulation scheme, we fix ni=4 measurements for each subject i and generate data according to the model of ( ² ) where β=(0.0,0.1,0.2,0.4,0.8) , bijk∼N(0,0.6) and εijk∼N(0,1) . 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.2,0.5,0.7,0.9 . Initial values of the variance components are set to be the estimates obtained from fitting a mixed model with no ridge component. In total, B=100 simulations are conducted for each of the ρ values based on n=40 individuals. For further assumed the expert prior knowledge is β=0 , for the computation of the PTMR and Stein-type estimators. Simulation results are summarized in Tables  ¹ and  ² . From Table  ² , it is clear the estimation is dependent on the level of significance α .

From Table  ¹ , 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  ² , as the level of significance increases, the MSE increases. The graphs of the MSE against the different values of ρ are also shown in Fig.  ¹ .

Fig. 1

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 [ ⁴ ], 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  ³ introduces variables of this study:

Variables have been measured ni 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. [ ¹³ ]). Information on the important variables is summarized in Table  ⁴ .

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  ⁵ shows the mixed, MR, shrinkage M , shrinkage MR and PTMR estimators, respectively, denoted by β^M , β^MR , β^MS , β^MRS and β^MRPT . From Table  ⁶ , 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. [ ²⁶ ]) 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 [ ²⁷ ] and Bilgrami and O’Keefe [ ²⁸ ].

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, {(Xtest,ytest)} , 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.,

Our results are based on N=500 case re-sampled bootstrap sample. In Table  ⁵ , 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.  ² .

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.