Sparse Estimation Strategies in Linear Mixed Effect Models for High-Dimensional Data Application

Abstract

In a host of business applications, biomedical and epidemiological studies, the problem of multicollinearity among predictor variables is a frequent issue in longitudinal data analysis for linear mixed models (LMM). We consider an efficient estimation strategy for high-dimensional data application, where the dimensions of the parameters are larger than the number of observations. In this paper, we are interested in estimating the fixed effects parameters of the LMM when it is assumed that some prior information is available in the form of linear restrictions on the parameters. We propose the pretest and shrinkage estimation strategies using the ridge full model as the base estimator. We establish the asymptotic distributional bias and risks of the suggested estimators and investigate their relative performance with respect to the ridge full model estimator. Furthermore, we compare the numerical performance of the LASSO-type estimators with the pretest and shrinkage ridge estimators. The methodology is investigated using simulation studies and then demonstrated on an application exploring how effective brain connectivity in the default mode network (DMN) may be related to genetics within the context of Alzheimer’s disease.

1. Introduction

In many fields such as bio-informatics, physical biology, and epidemiology, the response of interest is represented by repeated measures of some variables of interest that are collected over a specified time period for different independent subjects or individuals. These 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, and baseline measures, among others. A widely-used statistical tool in the analysis and modeling of longitudinal and repeated measures data is the linear mixed effects model (LMM) [1,2]. This model provides an effective and flexible way to describe the means and the covariance structures of a response variable after accounting for within subject correlation.

The rapid growth in the size and scope of longitudinal data has created a need for innovative statistical strategies in longitudinal data analysis. Classical methods are based on the assumption that the number of predictors is less than the number of observations. However, there is an increasing demand for efficient prediction strategies for analysis of high-dimensional data, where the number of observed data elements (sample size) are smaller than the number of predictors in a linear model context. Existing techniques that deal with high-dimensional data mostly rely on various penalized estimators. Due to the trade-off between model complexity and model prediction, the statistical inference of model selection becomes an extremely important and challenging problem in high-dimensional data analysis.

Over the years, many penalized regularization approaches have been developed to do variable selection and estimation simultaneously. Among them, the least absolute shrinkage and selection operator (LASSO) is commonly used [3]. It is a useful estimation technique in part due to its convexity and computational efficiency. The LASSO approach is based on an ${\ell }_1$ penalty for regularization of regression parameters. Ref. [4] provides a comprehensive summary of the consistency properties of the LASSO approach. Related penalized likelihood methods have been extensively studied in the literature, see for example [5,6,7,8,9,10]. The penalized likelihood methods have a close connection to Bayesian procedures. Thus, the LASSO estimate corresponds to a Bayes method that puts a Laplacian (double-exponential) prior on the regression coefficients [11,12].

In this paper, our interest lies in estimating the fixed effect parameters of the LMM using a ridge estimation technique when it is assumed that some prior information is available in the form of potential linear restrictions on the parameters. One possible source of prior information is using a Bayesian approach. An alternative source of prior information may be obtained from previous studies or expert knowledge that search for or assume sparsity patterns.

We consider the problem of fixed effect parameter estimation for LMMs when there exist many predictors relative to the sample size. These predictors may be classified into two groups: sparse and non-sparse. Thus, there are two choices to be considered: a full model with all predictors, and a sub-model that contains only non-sparse predictors. When the sub-model based on available subspace information is true (i.e., the assumed restriction holds), it then provides more efficient statistical inferences than those based on a full model. In contrast, if the sub-model is not true, the estimates could become biased and inefficient. The consequences of incorporating subspace information therefore depend on the quality or reliability of the information being incorporated into the estimation procedure. One way to deal with uncertain subspace information is to use a pretest estimation strategy. The validity of the information is tested before incorporation into a final estimator. Another approach is shrinkage estimation, which shrinks the full model estimator to the sub-model estimator by utilizing subspace information. Besides these estimation strategies, there is a growing literature on simultaneous model selection and estimation. These approaches are known as penalty strategies. By shrinking some regression coefficients toward zero, the penalty methods simultaneously select a sub-model and estimate its regression parameters. Several authors have investigated the pretest, shrinkage, and penalty estimation strategies in partial linear model, Poisson regression model, and Weibull censored regression model [13,14,15].

To formulate the problem, we suppose that the vector of the fixed effects parameter $\mathit{\beta }$ in the LMM can be partitioned into two sub-vectors $\mathit{\beta }={({\mathit{\beta }}_1^{\prime },{\mathit{\beta }}_2^{\prime })}^{\prime }$, where ${\mathit{\beta }}_1$ is the coefficient vector of non-sparse predictors and ${\mathit{\beta }}_2$ is the coefficient vector of sparse predictors. Our interest lies in the estimation of ${\mathit{\beta }}_1$ when ${\mathit{\beta }}_2$ is close to zero. To deal with this problem in the context of low dimensional data, ref. [16] propose an improved estimation strategy using sub-model selection and post-estimation for the LMM. Within this framework, linear shrinkage and shrinkage pretest estimation strategies are developed, which combine full model and sub-model estimators in an effective way as a trade-off between bias and variance. Ref. [17] extend this study by using a likelihood ratio test to develop James–Stein shrinkage and pretest estimation methods based on LMM for longitudinal data. In addition, the non-penalty estimators are compared with several penalty estimators (LASSO, adaptive LASSO and Elastic Net) for best performance.

In most real data situations, there is also the problem of multicollinearity among predictor variables for high-dimensional data. Various biased estimation techniques such as shrinkage estimation, partial least squares estimation [18] and Liu estimators [19] have been implemented to deal with this problem, but the widely used technique is ridge estimation [20]. The ridge estimator overcomes the weakness of the least squares estimator with a smaller mean squared error. To overcome and combat multicollinearity, ref. [21] propose pretest and Stein-type ridge regression estimators for linear and partially linear models. Furthermore, ref. [22] also develop shrinkage estimation based on Liu regression to overcome multicollinearity in linear models.

Our primary focus is on the estimation and prediction problem for linear mixed effect models when there are many potential predictors that have a weak or no influence on the response of interest. This method simultaneously controls overfitting using general least square estimation with a roughness penalty. We propose pretest and shrinkage estimation strategies using the ridge estimation technique as a base estimator and numerically compare their performance with the LASSO and adaptive LASSO estimators. Our proposed estimation strategy is applied to both high-dimensional and low-dimensional data.

The rest of this article is organized as follows. In Section 2, we present the linear mixed effect model and the proposed estimation techniques. We introduce the full and sub-model estimators based on ridge estimation. Thereafter, we construct the pretest and shrinkage ridge estimators. Section 3 provides the asymptotic bias and risk of these estimators. A Monte Carlo simulation is used to evaluate the performance of the estimators including a comparison with the lasso-type estimators, and the results are reported in Section 4. Section 5 presents a demonstration of the proposed methodology on a high-dimensional resting-state effective brain connectivity and genetic data. We also illustrate the proposed estimation methods in an application to a low-dimensional Amsterdam growth and health study. Section 6 presents a discussion with recommendations.

2. Model and Estimation Strategies

In this section, we present the linear mixed effect model and the proposed estimation strategies.

2.1. Linear Mixed Model

Suppose that we have a sample of N subjects. For the $i^{th}$ subject, we collect the response variable $y_{ij}$ for the jth time, where $i=1\dots ,n;j=1\dots ,n_i$ and $N={\sum }_{i=1}^n{n}_i$. Let ${\mathit{Y}}_i={(y_{i1},\dots y_{i{n}_i})}^{\prime }$ denotes the $n_i\times 1$ vector of responses from the ith subject. Let ${\mathbf{X}}_i={({\mathbf{x}}_{i1},\dots ,{\mathbf{x}}_{i{n}_i})}^{\prime }$ and ${\mathbf{Z}}_i={({\mathbf{z}}_{i1},\dots ,{\mathbf{z}}_{i{n}_i})}^{\prime }$ be $n_i\times$ p and $n_i\times$ q known fixed-effects and random-effect design matrix for the ith subject of full rank p and q, respectively. The linear mixed effect model [1] for a vector of repeated responses ${\mathit{Y}}_i$ on the ith subject is assumed to have the form

$${\mathbf{Y}}_i={\mathbf{X}}_i\mathit{\beta }+{\mathbf{Z}}_i{\mathbf{a}}_i+{\mathit{ϵ}}_i,$$

(1)

where $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_p)}^{\prime }$ is the p × 1 vector of unknown fixed-effect parameters or regression coefficients, ${\mathit{a}}_i$ is the q × 1 vector of unobservable random effects for the ith subject, assumed to come from a multivariate normal distribution with zero mean and a covariance matrix G, where G is an unknown $q\times q$ covariance matrix and ${ϵ}_i$ denotes $n_i\times$ 1 vector of error terms assumed to be normally distributed with zero mean, covariance matrix ${\sigma }^2{\mathbf{I}}_{n_i}$. Further, ${\mathit{ϵ}}_i$ are assumed to be independent of the random effects ${\mathit{a}}_i$.

The marginal distribution for the response $y_i$ is normal with mean ${\mathbf{X}}_i\mathit{\beta }$ and covariance matrix $Cov\left({\mathit{Y}}_i\right)={\mathit{Z}}_i{\sigma }_i^2{\mathit{Z}}_i^{\mathit{T}}+{\sigma }^2{\mathit{I}}_{\mathit{n}}.$ By stacking the vectors, the mixed model can be can be expressed as $\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathbf{Za}+\mathit{ϵ}$. From the Equation (1), the distribution of the model follows $\mathit{Y}\sim {\mathcal{N}}_n(\mathit{X}\mathit{\beta },\mathit{V})$, where $E\left(\mathit{Y}\right)=\mathit{X}\mathit{\beta }$ with covariance, $\mathit{V}=\sum _{i=1}^n{\mathit{Z}}_i{\sigma }_i^2{\mathit{Z}}_i^{\mathit{T}}+{\sigma }^2{\mathbf{I}}_n$.

2.2. Ridge Full Model and Sub-Model Estimator

The generalized least square estimator (GLS) is defined as ${\widehat{\mathit{\beta }}}^{\mathrm{GLS}}={\left({\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{Y}$ and the ridge full model estimator can be obtained by introducing a penalized regression so that $\widehat{\mathit{\beta }}=arg{min}_{\mathit{\beta }}\left\{{(\mathbf{Y}-\mathbf{X}\mathit{\beta })}^T{\mathbf{V}}^{-1}(\mathbf{Y}-\mathbf{X}\mathit{\beta })+k{\mathit{\beta }}^{\mathit{T}}\mathit{\beta }\right\}$ and ${\widehat{\mathit{\beta }}}^{\mathrm{Ridge}}={({\mathbf{X}}^T{\mathit{V}}^{-1}\mathbf{X}+k\mathbf{I})}^{-1}{\mathbf{X}}^T{\mathit{V}}^{-1}\mathbf{Y}$, where ${\widehat{\mathit{\beta }}}^{\mathrm{Ridge}}$ is the ridge full model estimator and $k\in [0,\infty )$ is the tuning parameter. If k = 0, ${\widehat{\mathit{\beta }}}^{\mathrm{Ridge}}$ is the GLS estimator and ${\widehat{\mathit{\beta }}}^{\mathrm{Ridge}}=0$ for k is sufficiently large. We select the value of k using cross validation.

We let $\mathbf{X}=({\mathbf{X}}_1,{\mathbf{X}}_2)$, where ${\mathbf{X}}_1$ is an $n\times p_1$ sub-matrix containing the non-sparse predictors and ${\mathbf{X}}_2$ is an $n\times p_2$ sub-matrix that contains the sparse predictors. Accordingly, $\mathit{\beta }=({\mathit{\beta }}_1,{\mathit{\beta }}_2)$ where ${\mathit{\beta }}_1$ and ${\mathit{\beta }}_2$ have dimensions $p_1$ and $p_2$, respectively, with $p_1+p_2=p$, $p_i\ge 0$ for i = 1, 2.

A sub-model is defined as $\mathit{Y}=\mathit{X}\mathit{\beta }+\mathit{Za}+\mathit{ϵ}\mathrm{subject}\mathrm{to}{\mathit{\beta }}^T\mathit{\beta }\le \varphi \mathrm{and}{\mathit{\beta }}_2=\mathbf{0}$ which corresponds to $\mathit{Y}={\mathit{X}}_{\mathbf{1}}{\mathit{\beta }}_{\mathbf{1}}+\mathit{Za}+\mathit{ϵ}\mathrm{subject}\mathrm{to}{\mathit{\beta }}_{\mathbf{1}}^T{\mathit{\beta }}_{\mathbf{1}}\le \varphi .$ The sub-model estimator ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}$ of ${\mathit{\beta }}_1$ has the form ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}={({\mathbf{X}}_1^T{\mathbf{V}}^{-1}{\mathbf{X}}_1+k\mathbf{I})}^{-1}{\mathbf{X}}_1^T{\mathbf{V}}^{-1}\mathbf{Y}.$ We denote ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$ as the full model ridge estimator of ${\mathit{\beta }}_1$ and given as ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}={({\mathbf{X}}_1^T{\mathbf{V}}^{-1/2}{\mathrm{M}}_{X_2}{\mathbf{V}}^{-1/2}{\mathbf{X}}_1+k\mathrm{I})}^{-1}{\mathbf{X}}_1^T{\mathbf{V}}^{-1/2}{\mathrm{M}}_{X_2}{\mathbf{V}}^{-1/2}\mathbf{Y},$ where ${\mathbf{M}}_{X_2}=\mathrm{I}-\mathbf{P}=\mathbf{I}-{\mathbf{V}}^{-1/2}{\mathbf{X}}_2{\left({\mathbf{X}}_2{\mathbf{V}}^{-1}{\mathbf{X}}_2\right)}^{-1}{\mathbf{X}}_2^T{\mathbf{V}}^{-1/2}.$

2.3. Pretest Ridge Estimation Strategy

Generally, the sub-model estimator will be more efficient than the full model estimator if the information embodied in the imposed linear restrictions is valid, thus ${\mathit{\beta }}_2$ is close to zero. However, if the information is not valid the sub-model estimator is likely to be more biased and may have a higher risk than the full model estimator. There is, therefore, some doubt as to whether or not to impose the restrictions on the model’s parameter. It is in response to this uncertainty that a statistical test may be used to determine the validity of the proposed restrictions. Accordingly, the procedure to follow in practice is pretest the validity of the restrictions and if the outcome of the pretest suggests that they are correct then the model parameters are estimated incorporating the restrictions. If the pretest rejects the restrictions then the parameters are estimated from the sample information alone. This motivates the consideration of the pretest estimation strategy for the LMM.

The pretest estimator is a combination of the full model estimator ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}},$ and sub-model estimator ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}},$ through an indicator function $\mathrm{I}({\mathrm{L}}_n\le d_{n,\alpha })$, where ${\mathrm{L}}_n$ is an appropriate test statistic to test $H_0:{\mathit{\beta }}_2=\mathbf{0}$ versus $H_A:{\mathit{\beta }}_2\ne \mathbf{0}$. Moreover, $d_{n,\alpha }$ is an $\alpha$ level critical value based on distribution of ${\mathrm{L}}_n$ under $H_0$. We define test statistics based on the log-likelihood ratio test as ${\mathrm{L}}_n=2\left\{{\ell }^{*}({\widehat{\mathit{\beta }}}^{\mathrm{RFM}}\mid \mathbf{Y})-{\ell }^{*}({\widehat{\mathit{\beta }}}^{\mathrm{RSM}}\mid \mathbf{Y})\right\}$.

Under $H_0$, the test statistic ${\mathrm{L}}_n$ follows asymptotic chi-square distribution with $p_2$ degrees of freedom. The pretest test ridge estimator ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}}$ of ${\mathit{\beta }}_1$ is then defined by

$${\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}}={\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}-({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}-{\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}})\mathrm{I}({\mathrm{L}}_n\le d_{n,\alpha }),p_2\ge 1$$

.

2.4. Shrinkage Ridge Estimation Strategy

The pre-test estimator is a discontinuous function of the sub-model ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}$ and full model ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$, which depends on the hard threshold $(d_{n,\alpha }={\chi }_{p_2,\alpha }^2)$. We address this limitation by defining the shrinkage ridge estimator based on soft thresholding. The shrinkage ridge estimator (RSE) of ${\mathit{\beta }}_1$, denoted as ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}},$ is defined as

$${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}}={\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}+({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}-{\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}})(1-(p_2-2){\mathrm{L}}_n^{-1}),p_2\ge 3.$$

Here, ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}}$ is the linear combination of the full model ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$ and sub-model ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}$ estimates. If ${\mathrm{L}}_n\le (p_2-2),$ then a relatively large weight is placed on ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}$ otherwise, more weight is on ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$. A setback with ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}}$ is that it is not a convex combination of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$ and ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}$. This can cause over-shrinkage, which gives the estimator opposite sign of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$. This could happen if $(p_2-2){\mathrm{L}}_n^{-1}$ is larger than one. To counter this, we use the positive-part shrinkage ridge estimator (RPS) defined as

$${\widehat{\mathit{\beta }}}_1^{\mathrm{RPS}}={\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}+({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}-{\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}){(1-(p_2-2)L_n^{-1})}^{+},p_2\ge 3$$

where ${(1-(p_2-2){\mathrm{L}}_n^{-1})}^{+}=max(0,1-(p_2-2){\mathrm{L}}_n^{-1}).$ The RPS estimator will control possible over-shrinking in the RSE estimator.

3. Asymptotic Results

In this section, we derive the asymptotic distributional bias and risk of the estimators considered in Section 2. We examine the properties of the estimators for increasing n and as ${\mathit{\beta }}_2$ approaches the null vector under the sequence of local alternatives defined as

$$K_n:{\mathit{\beta }}_2={\mathit{\beta }}_{2\left(n\right)}=\frac{\mathit{\kappa }}{\sqrt{n}},$$

(2)

where $\mathit{\kappa }={({\kappa }_1,{\kappa }_2\dots ,{\kappa }_{p_2})}^{\prime }\in {\mathbb{R}}^{p_2}$ is a fixed vector. The vector $\frac{\mathit{\kappa }}{\sqrt{n}}$ is a measure of how far local alternatives $K_n$ differ from the subspace information ${\mathit{\beta }}_2=\mathbf{0}$. In order to evaluate the performance of the estimators, we define the asymptotic distributional bias of the estimator ${\widehat{\mathit{\beta }}}_1^{*}$ as

$$\mathrm{ADB}\left({\widehat{\mathit{\beta }}}_1^{*}\right)=\underset{n\to \infty }{lim}E\left\{\sqrt{n}({\widehat{\mathit{\beta }}}_1^{*}-{\mathit{\beta }}_1)\right\},$$

In order to compute the risk functions, we first compute the asymptotic covariance of the estimators. The asymptotic covariance of an estimator ${\widehat{\mathit{\beta }}}_1^{*}$ is expressed as

$$\mathrm{Cov}\left({\widehat{\mathit{\beta }}}_1^{*}\right)=\underset{n\to \infty }{lim}E\left\{n({\widehat{\mathit{\beta }}}_1^{*}-{\mathit{\beta }}_1){({\widehat{\mathit{\beta }}}_1^{*}-{\mathit{\beta }}_1)}^{\mathrm{T}}\right\}.$$

Following the asymptotic covariance matrix, we define the asymptotic risk of an estimator ${\widehat{\mathit{\beta }}}_1^{*}$ as $\mathrm{R}\left({\widehat{\mathit{\beta }}}_{\mathbf{1}}^{*}\right)=\mathrm{tr}\left(\mathbf{Q}\mathrm{Cov}\left({\widehat{\mathit{\beta }}}_1^{*}\right)\right)$. Q is a positive definite matrix of weights with dimensions of $p\times p$. We set Q = I in this study.

Assumption 1.

We make the following two regularity conditions to establish the asymptotic properties of the estimators.

1.$\frac{1}{n}{max}_{1\le i\le n}^{}{\mathit{x}}_i^T{\left[{\mathit{X}}^T{\mathit{V}}^{-1}\mathit{X}\right]}^{-1}{\mathit{x}}_i\to$0as n$\to \infty ,$where${\mathit{x}}_i^T$is the ith row ofX.

2.${\mathit{B}}_n=n^{-1}{\left[{\mathit{X}}^T{\mathit{V}}^{-1}\mathit{X}\right]}^{-1}\to \mathit{B},$for some finiteB = $\left(\begin{array}{cc}{\mathit{B}}_{11}& {\mathit{B}}_{12}\\ {\mathit{B}}_{21}& {\mathit{B}}_{22}\end{array}\right)$.

Theorem 1.

For$k<\infty$, If$k/\sqrt{n}\to {\lambda }_o$and B is non-singular, the distribution of the full model ridge estimator, ${\widehat{\mathit{\beta }}}_n^{\mathrm{RFM}}$is

$$\sqrt{n}({\widehat{\mathit{\beta }}}_n^{\mathrm{RFM}}-\mathit{\beta })\stackrel{D}{\to }\mathcal{N}(-{\lambda }_o{\mathit{B}}^{-1}\mathit{\beta },{\mathit{B}}^{-1}),$$

where$\stackrel{D}{\to }$denotes convergence in distribution.

Proof. 

See Theorem 2 in [23]. □

Proposition 1.

Assuming the above assumption 1 together with Theorem 1 hold, under the local alternatives$K_n$, we have

$$\left(\begin{array}{c}{\phi }_1\\ {\phi }_3\end{array}\right)\stackrel{D}{\to }\mathcal{N}\left[\left(\begin{array}{c}-{\mathit{\mu }}_{11.2}\\ \mathit{\delta }\end{array}\right),\left(\begin{array}{cc}{\mathit{B}}_{11.2}^{-1}& \mathbf{\Phi }\\ \mathbf{\Phi }& \mathbf{\Phi }\end{array}\right)\right],$$

$$\left(\begin{array}{c}{\phi }_3\\ {\phi }_2\end{array}\right)\stackrel{D}{\to }\mathcal{N}\left[\left(\begin{array}{c}\mathit{\delta }\\ -\gamma \end{array}\right),\left(\begin{array}{cc}\mathbf{\Phi }& \mathbf{0}\\ \mathbf{0}& {\mathbf{B}}_{11}^{-1}\end{array}\right)\right],$$

where ${\phi }_1=\sqrt{n}({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}-{\mathit{\beta }}_1)$, ${\phi }_2=\sqrt{n}({\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}-{\mathit{\beta }}_1)$, ${\phi }_3=\sqrt{n}({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}-{\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}})$, $\gamma ={\mathit{\mu }}_{11.2}+\mathit{\delta }$, $\mathit{\delta }={\mathbf{B}}_{11}^{-1}{\mathbf{B}}_{12}\mathit{\kappa }$, $\mathbf{\Phi }={\mathbf{B}}_{11}^{-1}{\mathbf{B}}_{12}{\mathbf{B}}_{22.1}^{-1}{\mathrm{B}}_{21}{\mathbf{B}}_{11}^{-1}$, ${\mathbf{B}}_{22.1}={\mathbf{B}}_{22}-{\mathbf{B}}_{21}{\mathbf{B}}_{11}^{-1}{\mathbf{B}}_{12}$, $\mathit{\mu }=-{\lambda }_o{\mathbf{B}}^{-1}\mathit{\beta }=\left(\begin{array}{c}{\mathit{\mu }}_1\\ {\mathit{\mu }}_2\end{array}\right)$ and ${\mathit{\mu }}_{11.2}={\mathit{\mu }}_1-{\mathbf{B}}_{12}{\mathbf{B}}_{22}^{-1}(({\mathit{\beta }}_2-\mathit{\kappa })-{\mathit{\mu }}_2)$.

Proof. 

See Appendix A. □

Theorem 2.

Under the condition of Theorem 1 and the local alternatives$K_n$, the ADBs of the proposed estimators are

$$\begin{array}{cc}\hfill ADB\left({\widehat{\mathit{\beta }}}_1^{\mathit{RFM}}\right)& =-{\mathit{\mu }}_{11.2},\hfill \\ \hfill ADB\left({\widehat{\mathit{\beta }}}_1^{\mathit{RSM}}\right)& =-{\mathit{\mu }}_{11.2}-{\mathit{B}}_{11}^{-1}{\mathit{B}}_{12}\mathit{\delta }=-\gamma ,\hfill \\ \hfill ADB\left({\widehat{\mathit{\beta }}}_1^{\mathit{RPT}}\right)& =-{\mathit{\mu }}_{11.2}-\mathit{\delta }{\mathit{H}}_{p_2+2}({\chi }_{p_2,\alpha }^2;\Delta ),\hfill \\ \hfill ADB\left({\widehat{\mathit{\beta }}}_1^{\mathit{RSE}}\right)& =-{\mathit{\mu }}_{11.2}-(p_2-2)\mathit{\delta }E\left({\chi }_{p_2+2}^{-2}(\Delta )\right),\hfill \\ \hfill ADB\left({\widehat{\mathit{\beta }}}_1^{\mathit{RPS}}\right)& =-{\mathit{\mu }}_{11.2}-\mathit{\delta }{\mathit{H}}_{p_2+2}({\chi }_{p_2-2}^2;\Delta )\}-(p_2-2)\mathit{\delta }E\left\{{\chi }_{p_2+2}^{-2}(\Delta )I({\chi }_{p_2+2}^{-2}>p_2-2)\right\},\hfill \end{array}$$

where$\Delta ={\mathit{\kappa }}^T{\mathit{B}}_{22.1}^{-1}\mathit{\kappa }$, ${\mathit{B}}_{22.1}={\mathit{B}}_{22}-{\mathit{B}}_{21}{\mathit{B}}_{11}^{-1}{\mathit{B}}_{12}$, and $H_v(x;\Delta )$is the cumulative distribution function of the non-central chi-squared distribution with non-centrality parameter Δ and v degrees of freedom, and $E\left({\chi }_v^{-2j}(\Delta )\right)$ is the expected value of the inverse of a non-central ${\chi }^2$ distribution with v degrees of freedom and non-centrality parameter Δ,

$$E\left({\chi }_v^{-2j}(\Delta )\right)={\int }_0^{\infty }{x}^{-2j}d{H}_v(x,\Delta ).$$

Proof. 

See Appendix B.1. □

Since the ADBs of the estimators are in non-scalar form, we define the following asymptotic quadratic bias (AQDB) of ${\widehat{\mathit{\beta }}}_1^{*}$ by

$$\mathrm{AQDB}\left({\widehat{\mathit{\beta }}}_1^{*}\right)={\left(\mathrm{ADB}\left({\widehat{\mathit{\beta }}}_1^{*}\right)\right)}^{\prime }{\mathbf{B}}_{11.2}\left(\mathrm{ADB}\left({\widehat{\mathit{\beta }}}_1^{*}\right)\right),$$

where ${\mathbf{B}}_{11.2}={\mathbf{B}}_{11}-{\mathbf{B}}_{12}{\mathbf{B}}_{22}^{-1}{\mathbf{B}}_{21}.$

Corollary 1.

Suppose Theorem 2 holds. Then, under$\left\{K_n\right\}$, the AQDBs of the estimators are

$$\begin{array}{cc}\hfill \mathit{AQDB}\left({\widehat{\mathit{\beta }}}_1^{\mathit{RFM}}\right)& ={\mathit{\mu }}_{11.2}^T{\mathit{B}}_{11.2}{\mathit{\mu }}_{11.2},\hfill \\ \hfill \mathit{AQDB}\left({\widehat{\mathit{\beta }}}_1^{\mathit{RSM}}\right)& ={\gamma }^T{\mathit{B}}_{11.2}\gamma ,\hfill \\ \hfill \mathit{AQDB}\left({\widehat{\mathit{\beta }}}_1^{\mathit{RPT}}\right)& ={\mathit{\mu }}_{11.2}^T{\mathit{B}}_{11.2}{\mathit{\mu }}_{11.2}+{\mathit{\mu }}_{11.2}^T{\mathit{B}}_{11.2}\mathit{\delta }{\mathit{H}}_{p_2+2}({\chi }_{p_2}^2;\Delta )\hfill \\ & +{\mathit{\delta }}^T{\mathit{B}}_{11.2}{\mathit{\mu }}_{11.2}{\mathit{H}}_{p_2+2}({\chi }_{p_2}^2;\Delta )+{\mathit{\delta }}^T{\mathit{B}}_{11.2}\mathit{\delta }{\mathit{H}}_{p_2+2}^2({\chi }_{p_2}^2;\Delta ),\hfill \\ \hfill \mathit{AQDB}\left({\widehat{\mathit{\beta }}}_1^{\mathit{RSE}}\right)& ={\mathit{\mu }}_{11.2}^T{\mathit{B}}_{11.2}{\mathit{\mu }}_{11.2}+(p_2-2){\mathit{\mu }}_{11.2}^T{\mathit{B}}_{11.2}\mathit{\delta }E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)\hfill \\ & +(p_2-2){\mathit{\delta }}^T{\mathit{B}}_{11.2}{\mathit{\mu }}_{11.2}E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)+{(p_2-2)}^2{\mathit{\delta }}^T{\mathit{B}}_{11.2}\mathit{\delta }{\left(E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)\right)}^2,\hfill \\ \hfill \mathit{AQDB}\left({\widehat{\mathit{\beta }}}_1^{\mathit{RPS}}\right)& ={\mathit{\mu }}_{11.2}^T{\mathit{B}}_{11.2}{\mathit{\mu }}_{11.2}+\left({\mathit{\delta }}^T{\mathit{B}}_{11.2}{\mathit{\mu }}_{11.2}+{\mathit{\mu }}_{11.2}^T{\mathit{B}}_{11.2}\mathit{\delta }\right)[{\mathit{H}}_{p_2+2}(p_2-2;\Delta )\hfill \\ & +(p_2-2)E\left\{{\chi }_{p_2+2}^{-2}(\Delta )I({\chi }_{p_2+2}^{-2}(\Delta )>p_2-2)\right\}]+{\mathit{\delta }}^T{\mathit{B}}_{11.2}\mathit{\delta }[{\mathit{H}}_{p_2+2}(p_2-2;\Delta )\hfill \\ & +(p_2-2)E\left\{{\chi }_{p_2+2}^{-2}(\Delta )I({\chi }_{p_2+2}^{-2}(\Delta )>p_2-2)\right\}{]}^2.\hfill \end{array}$$

When ${\mathbf{B}}_{11.2}=\mathbf{0}$, the AQDB of all estimators are equivalent, and the estimators are therefore asymptotically unbiased. If we assume that ${\mathbf{B}}_{11.2}\ne 0$, the results for the bias of the estimators can be summarized as follows:

• The AQDB of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}$ is an unbounded function of ${\gamma }^T{\mathbf{B}}_{11.2}\gamma$.
• The AQDB of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}}$ starts from ${\mathit{\mu }}_{11.2}^T{\mathbf{B}}_{11.2}{\mathit{\mu }}_{11.2}$ at $\Delta =0$, and when $\Delta$ increases, it increases to the maximum and then decreases to zero.
• The characteristics of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}}$ and ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPS}}$ are similar to ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}}$. The AQDB of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}}$ and ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPS}}$ similarly start from ${\mathit{\mu }}_{11.2}^T{\mathbf{B}}_{11.2}{\mathit{\mu }}_{11.2}$ at $\Delta =0$, and increase to a point, and then decrease towards zero, since $E\left\{{\chi }_{p_2+2}^{-2}(\Delta )\right\}$ is a non-increasing on of $\Delta$.

Theorem 3.

Suppose Theorem 1 holds and under the local alternatives$K_n$, the covariance matrices of the estimators are

$$\begin{array}{cc}\hfill Cov\left({\widehat{\mathit{\beta }}}_{\mathbf{1}}^{\mathit{RFM}}\right)& ={\mathit{B}}_{11.2}^{-1}+{\mathit{\mu }}_{11.2}{\mathit{\mu }}_{11.2}^T,\hfill \\ \hfill Cov\left({\widehat{\mathit{\beta }}}_{\mathbf{1}}^{\mathit{RSM}}\right)& ={\mathit{B}}_{11}^{-1}+\gamma {\gamma }^T,\hfill \\ \hfill Cov\left({\widehat{\mathit{\beta }}}_{\mathbf{1}}^{\mathit{RPT}}\right)& ={\mathit{B}}_{11.2}^{-1}+{\mathit{\mu }}_{11.2}{\mathit{\mu }}_{11.2}^T+2{\mathit{\mu }}_{11.2}^T\mathit{\delta }{\mathit{H}}_{p_2+2}({\chi }_{p_2}^2;\Delta )-\mathbf{\Phi }{\mathit{H}}_{p_2+2}({\chi }_{p_2}^2;\Delta )\hfill \\ & +\mathit{\delta }{\mathit{\delta }}^{\mathit{T}}\left[2{\mathit{H}}_{p_2+2}({\chi }_{p_2}^2;\Delta )-{\mathit{H}}_{p_2+4}({\chi }_{p_2}^2;\Delta )\right],\hfill \\ \hfill Cov\left({\widehat{\mathit{\beta }}}_{\mathbf{1}}^{\mathit{RSE}}\right)& ={\mathit{B}}_{11.2}^{-1}+{\mathit{\mu }}_{11.2}{\mathit{\mu }}_{11.2}^T+2(p_2-2){\mathit{\mu }}_{11.2}^T\mathit{\delta }E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)\hfill \\ & -(p_2-2)\mathbf{\Phi }\left\{2E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)-(p_2-2)E\left({\chi }_{p_2+2}^{-4}(\Delta )\right)\right\}\hfill \\ & +(p_2-2)\mathit{\delta }{\mathit{\delta }}^T\left\{-2E\left({\chi }_{p_2+4}^{-2}(\Delta )\right)+2E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)+(p_2-2)E\left({\chi }_{p_2+4}^{-4}(\Delta )\right)\right\},\hfill \\ \hfill Cov\left({\widehat{\mathit{\beta }}}_1^{\mathit{RPS}}\right)& =\mathit{Cov}\left({\widehat{\mathit{\beta }}}_1^{\mathit{RSE}}\right)+2\mathit{\delta }{\mathit{\mu }}_{11.2}^{T}E\left(\left\{1-(p_2-2){\chi }_{p_2+2}^{-2}(\Delta )\right\}I\left({\chi }_{p_2+2}^2(\Delta )\le p_2-2\right)\right)\hfill \\ & -2\mathbf{\Phi }E\left(\left\{1-(p_2-2){\chi }_{p_2+2}^{-2}(\Delta )\right\}I\left({\chi }_{p_2+2}^2(\Delta )\le p_2-2\right)\right)\hfill \\ & -2\mathit{\delta }{\mathit{\delta }}^{T}E\left(\{1-(p_2-2){\chi }_{p_2+4}^{-2}(\Delta )\}I({\chi }_{p_2+4}^2(\Delta )\le p_2-2)\right)\hfill \\ & +2\mathit{\delta }{\mathit{\delta }}^{T}E\left(\left\{1-(p_2-2){\chi }_{p_2+2}^{-2}(\Delta )\right\}I\left({\chi }_{p_2+2}^2(\Delta )\le p_2-2\right)\right)\hfill \\ & -{(p_2-2)}^2\mathbf{\Phi }E\left({\chi }_{p_2+2}^{-4}(\Delta )I\left({\chi }_{p_2+2,\alpha }^2(\Delta )\le p_2-2\right)\right)\hfill \\ & -{(p_2-2)}^2\mathit{\delta }{\mathit{\delta }}^{T}E\left({\chi }_{p_2+2,\alpha }^{-4}(\Delta )I\left({\chi }_{p_2+2,\alpha }^2(\Delta )\le p_2-2\right)\right)\hfill \\ & +\mathbf{\Phi }{\mathit{H}}_{p_2+2}\left(p_2-2;\Delta \right)+\mathit{\delta }{\mathit{\delta }}^T{\mathit{H}}_{p_2+4}\left(p_2-2;\Delta \right).\hfill \end{array}$$

Proof. 

See Appendix B.2. □

Corollary 2.

Under the local alternatives ($K_n$) and from Theorem 3, the risk of the estimators are obtained as

$$\begin{array}{cc}\hfill R\left[{\widehat{\mathit{\beta }}}_1^{\mathit{RFM}}\right]& =tr\left({\mathit{Q}\mathit{B}}_{11.2}^{-1}\right)+{\mathit{\mu }}_{11.2}^T\mathit{Q}{\mathit{\mu }}_{11.2},\hfill \\ \hfill R[{\widehat{\mathit{\beta }}}_1^{\mathit{RSM}}]& =tr\left({\mathit{Q}\mathit{B}}_{11}^{-1}\right)+{\gamma }^T\mathit{Q}\gamma ,\hfill \\ \hfill R\left[{\widehat{\mathit{\beta }}}_1^{\mathit{RPT}}\right]& =tr\left({\mathit{Q}\mathit{B}}_{11.2}^{-1}\right)+{\mathit{\mu }}_{11.2}^T\mathit{Q}{\mathit{\mu }}_{11.2}+2{\mathit{\mu }}_{11.2}^T\mathit{Q}\mathit{\delta }{\mathit{H}}_{p_2+2}\left({\chi }_{p_2}^2;\Delta \right)\hfill \\ & -tr\left(\mathit{Q}\mathbf{\Phi }\right){\mathit{H}}_{p_2+2}\left({\chi }_{p_2}^2;\Delta \right)+\mathit{\delta }Q{\mathit{\delta }}^T\left[2{\mathit{H}}_{p_2+2}\left({\chi }_{p_2}^2;\Delta \right)-{\mathit{H}}_{p_2+4}\left({\chi }_{p_2}^2;\Delta \right)\right],\hfill \\ \hfill R\left[{\widehat{\mathit{\beta }}}_1^{\mathit{RSE}}\right]& =tr\left({\mathit{Q}\mathit{B}}_{11.2}^{-1}\right)+{\mathit{\mu }}_{11.2}^T\mathit{Q}{\mathit{\mu }}_{11.2}+2(p_2-2){\mathit{\mu }}_{11.2}^T\mathit{Q}\mathit{\delta }E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)\hfill \\ & -(p_2-2)\mathit{tr}\left(\mathit{Q}\mathbf{\Phi }\right)\left[E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)-(p_2-2)E\left({\chi }_{p_2+2}^{-4}(\Delta )\right)\right]\hfill \\ & +(p_2-2){\mathit{\delta }}^T\mathit{Q}\mathit{\delta }\left[2E\left({\chi }_{p_2+2}^{-2}(\Delta )\right)-2E\left({\chi }_{p_2+4}^{-2}(\Delta )\right)-(p_2-2)E\left({\chi }_{p_2+4}^{-4}(\Delta )\right)\right],\hfill \\ \hfill R\left[{\widehat{\mathit{\beta }}}_1^{\mathit{RPS}}\right]& =R\left[{\widehat{\mathit{\beta }}}_1^{\mathit{RSE}}\right]+2\mathit{\delta }\mathit{Q}{\mathit{\mu }}_{11.2}^{T}E\left(\left\{1-(p_2-2){\chi }_{p_2+2}^{-2}(\Delta )\right\}\mathrm{I}\left({\chi }_{p_2+2}^2(\Delta )\le p_2-2\right)\right)\hfill \\ & -2tr\left(\mathit{Q}\mathbf{\Phi }\right)E\left(\left\{1-(p_2-2){\chi }_{p_2+2}^{-2}(\Delta )\right\}I\left({\chi }_{p_2+2}^2(\Delta )\le p_2-2\right)\right)\hfill \\ & -2{\mathit{\delta }}^T\mathit{Q}\mathit{\delta }E\left(\{1-(p_2-2){\chi }_{p_2+4}^{-2}(\Delta )\}I({\chi }_{p_2+4}^2(\Delta )\le p_2-2)\right)\hfill \\ & +2{\mathit{\delta }}^T\mathit{Q}\mathit{\delta }E\left(\left\{1-(p_2-2){\chi }_{p_2+2}^{-2}(\Delta )\right\}I\left({\chi }_{p_2+2}^2(\Delta )\le p_2-2\right)\right)\hfill \\ & -{(p_2-2)}^{2}tr\left(\mathit{Q}\mathbf{\Phi }\right)E\left({\chi }_{p_2+2}^{-4}(\Delta )I\left({\chi }_{p_2+2}^2(\Delta )\le p_2-2\right)\right)\hfill \\ & -{(p_2-2)}^2{\mathit{\delta }}^T\mathit{Q}\mathit{\delta }E\left({\chi }_{p_2+2}^{-4}(\Delta )I\left({\chi }_{p_2+2}^2(\Delta )\le p_2-2\right)\right)\hfill \\ & +tr\left(\mathit{Q}\mathbf{\Phi }\right){\mathit{H}}_{p_2+2}\left(p_2-2;\Delta \right)+{\mathit{\delta }}^T\mathit{Q}\mathit{\delta }{\mathit{H}}_{p_2+4}\left(p_2-2;\Delta \right).\hfill \end{array}$$

From Theorem 2, when ${\mathbf{B}}_{12}=\mathbf{0},$ the risks of estimators ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}},$${\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}},$${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}},$ and ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPS}}$ are reduced to common value $\mathrm{tr}\left({\mathbf{Q}\mathbf{B}}_{11.2}^{-1}\right)+{\mathit{\mu }}_{11.2}^T\mathbf{Q}{\mathit{\mu }}_{11.2}$, the risk of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$. If ${\mathbf{B}}_{12}\ne \mathbf{0}$, the results can be summarized as follows:

• The risk of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$ remains constant while the risk of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}$ is an unbounded function of $\Delta$ since $\Delta \in [0,\infty ).$
• The risk of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}}$ increases as $\Delta$ moves away from zero, achieves it maximum and then decreases towards the risk of the full model estimator.
• The risk of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$ is smaller than the risk of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}}$ for small values in the neighborhood of $\Delta$ and for the rest of the parameter space, ${\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}}$ outperforms ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$, thus, $\mathrm{R}\left[{\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}\right]>\mathrm{R}\left[{\widehat{\mathit{\beta }}}_1^{\mathrm{RPT}}\right]$.
• Comparing the risks of ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}}$ and ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$, it can be seen that the estimator ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}}$ outperforms ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$ that is, $\mathrm{R}\left[{\widehat{\mathit{\beta }}}_1^{\mathrm{RSE}}\right]\le \mathrm{R}\left[{\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}\right]$ for all $\Delta \ge 0.$

4. Simulation Studies

In this section, we conduct a simulation study to assess the performance of the suggested estimators for finite samples. The criterion for comparing the performance of any estimator in our study is the mean square error. We simulate the response from the following LMM model

$${\mathit{Y}}_i={\mathbf{X}}_i\mathit{\beta }+{\mathbf{Z}}_i{\mathbf{a}}_i+{\mathit{ϵ}}_i,$$

(3)

where ${\mathit{ϵ}}_i\sim \mathcal{N}(0,{\sigma }^2{\mathrm{I}}_{n_i})$ with ${\sigma }^2=1$. We generate random effect covariate ${\mathbf{a}}_i$ from a multivariate normal distribution with zero mean and covariance matrix $\mathbf{G}=0.5{\mathrm{I}}_{2\times 2}$, where ${\mathrm{I}}_{2\times 2}$ is $2\times 2$ identity matrix. The design matrix ${\mathbf{X}}_i={({\mathbf{x}}_{i1},\dots ,{\mathbf{x}}_{i{n}_i})}^{\prime }$ is generated from a $n_i$-multivariate normal distribution with mean vector and covariance matrix ${\sum }_x$. Furthermore, we assume that the off-diagonal elements of the covariance matrix ${\sum }_x$ are equal to $\rho$, which is the coefficient of correlation between any two predictors, with $\rho =0.3,0.7,0.9.$ The ratio of the largest eigenvalue to the smallest eigen-value of matrix ${\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{X}$ is calculated as a condition number index (CNI) [24], which assesses the existence of multicollinearity in the design matrix. If the CNI is larger than 30, then the model has significant multicollinearity. Our simulations are based on the linear mixed effects model in Equation (3) with $n=60$ and 100 subjects.

We consider a situation when the model is assumed to be sparse. In this study, our interest lies in testing the hypothesis $H_o:{\mathit{\beta }}_2=\mathbf{0}$, and our goal is to estimate the fixed effect coefficient ${\mathit{\beta }}_1$. We partition the fixed effects coefficients as $\mathit{\beta }={({\mathit{\beta }}_1^{\prime },{\mathit{\beta }}_2^{\prime })}^{\prime }={({\mathit{\beta }}_1^{\prime },{\mathbf{0}}_{p_2})}^{\prime }.$ The coefficients ${\mathit{\beta }}_1$ and ${\mathit{\beta }}_2$ are $p_1$ and $p_2$ dimensional vectors, respectively, with $p=p_1+p_2$.

In order to investigate the behavior of the estimators, we define ${\Delta }^{*}=\left|\right|\mathit{\beta }-{\mathit{\beta }}_o\left|\right|$, where ${\mathit{\beta }}_o={({\beta }_1^T,{\mathbf{0}}_{p_2})}^{\mathrm{T}}$ and $\left|\right|.\left|\right|$ is the euclidean norm. We considered ${\Delta }^{*}$ values between 0 and 4. If ${\Delta }^{*}=0$, then we will have $\mathit{\beta }={(1,1,1,1,\underset{p_2}{\underbrace{0,0,\dots ,0}})}^{\mathrm{T}}$ to generate the response under null hypothesis. On the other hand, when ${\Delta }^{*}\ge 0,$ say ${\Delta }^{*}=4,$ we will have $\mathit{\beta }={(1,1,1,1,4,\underset{p_2-1}{\underbrace{0,0,\dots ,0}})}^{\mathrm{T}}$ to generate the response under the local alternative hypothesis. In our simulation study, we consider the number of fixed effect or predictor variables as $(p_1,p_2)\in \{(5,40),(5,500),(5,1000)\}.$ Each realization is repeated 5000 times to obtain consistent results and compute the MSE of suggested estimators with $\alpha =0.05$.

Based on the simulated data, we calculate the mean square error (MSE) of all the estimators as MSE$\left(\widehat{\mathit{\beta }}\right)=\frac{1}{5000}{\sum }_{j=1}^{5000}{(\widehat{\beta }-\beta )}^{\mathrm{T}}(\widehat{\beta }-\beta )$, where $\widehat{\mathit{\beta }}$ denotes any one of ${\widehat{\mathit{\beta }}}^{\mathrm{RSM}},{\widehat{\mathit{\beta }}}^{\mathrm{RPT}},{\widehat{\mathit{\beta }}}^{\mathrm{RSE}}$ and ${\widehat{\mathit{\beta }}}^{\mathrm{RPS}}$, in the jth repetition. We use the relative mean squared efficiency (RMSE), or the ratio of MSE for risk performance comparison. The RMSE of an estimator ${\widehat{\mathit{\beta }}}^{*}$ with respect to the baseline full model ridge estimator ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}$ is defined as $\mathrm{RMSE}({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}:{\widehat{\mathit{\beta }}}_1^{*})={\displaystyle \frac{\mathrm{MSE}\left({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}\right)}{\mathrm{MSE}\left({\widehat{\mathit{\beta }}}_1^{*}\right)}},$ where ${\mathit{\beta }}_1^{*}$ is one of the suggested estimators under consideration.

4.1. Simulation Results

In this subsection, we present the results from our simulation study. We report the results for $n=60,100$ and $p_1=5$ with different values of correlation coefficient $\rho$ are shown in

Table 1. RMSEs of RSM, RPT, RSE, and RPS estimators with respect to ${\widehat{\mathit{\beta }}}^{\mathrm{RFM}}$ when $\Delta \ge 0$ for $p_1=5$ and $n=60$.

$\mathit{\rho }$	${\mathit{p}}_2$	$\mathbf{\Delta }$	CNI	RSM	RPT	RSE	RPS
0.3	40	0	361	2.61	2.07	1.94	1.96
		1		1.05	1.07	1.20	1.25
		2		0.25	0.95	1.04	1.05
		3		0.12	0.98	0.99	1.00
		4		0.08	1.00	1.00	1.00
	500	0	613	4.48	3.29	3.48	1.96
		1		1.26	1.12	1.26	1.29
		2		0.41	0.97	1.08	1.09
		3		0.18	0.99	1.00	1.00
		4		0.13	1.00	1.00	1.00
	1000	0	693	5.36	4.53	4.67	4.71
		1		1.53	1.21	1.35	1.39
		2		0.49	1.01	1.13	1.14
		3		0.28	0.99	0.99	0.99
		4		0.10	1.00	1.00	1.00
0.7	40	0	1352	3.18	2.33	2.17	2.18
		1		1.04	1.11	1.20	1.23
		2		0.42	1.03	1.04	1.04
		3		0.23	0.98	0.99	1.00
		4		0.14	1.00	1.00	1.00
	500	0	1789	4.48	2.76	2.94	3.02
		1		1.08	1.43	1.52	1.53
		2		0.67	1.03	1.07	1.06
		3		0.35	0.98	1.00	1.00
		4		0.19	1.00	1.00	1.00
	1000	0	2134	6.82	5.24	5.30	3.02
		1		1.16	1.32	1.42	1.53
		2		0.75	1.10	1.15	1.16
		3		0.39	0.99	1.00	1.00
		4		0.11	1.00	1.00	1.00

. Furthermore, we plot the RMSEs against ${\Delta }^{*}$ in Figure 1 and Figure 2. The findings can be summarized as follows:

• When ${\Delta }^{*}=0$, the sub-model RSM outperforms all other estimators. As ${\Delta }^{*}=0$ moves from zero, the RMSE of the sub-model decreases and goes to zero.
• The pretest ridge estimator RPT outperforms shrinkage ridge and positive Stein ridge estimators in the case of ${\Delta }^{*}=0$. However, for large number of sparse predictors $p_2$ while keeping $p_1$ and n fixed, RPT is less efficient than RPS and RSE. In the case of ${\Delta }^{*}$ being larger than zero, the RMSE of RPT decreases, and it remains below 1 for immediate values of ${\Delta }^{*}$, after that the RMSE of RPT increases and approaches one for larger values of ${\Delta }^{*}.$
• RPS performs better than RSE in the entire parameter space induced by ${\Delta }^{*}$ as presented in Table 1 and

  Table 2. RMSEs of RSM, RPT, RSE, and RPS estimators with respect to ${\widehat{\mathit{\beta }}}^{\mathrm{RFM}}$ when $\Delta \ge 0$ for $p_1=5$, and $n=100$.

  $\mathit{\rho }$	${\mathit{p}}_2$	$\mathbf{\Delta }$	CNI	RSM	RPT	RSE	RPS
  0.3	40	0	150	2.38	2.09	1.88	1.90
  		1		0.89	1.01	1.05	1.08
  		2		0.21	0.94	1.01	1.02
  		3		0.06	0.94	0.99	1.00
  		4		0.02	1.00	1.00	1.00
  	500	0	340	4.15	2.65	2.99	3.17
  		1		0.87	1.08	1.18	1.21
  		2		0.14	0.96	1.03	1.05
  		3		0.06	0.99	0.99	1.00
  		4		0.03	1.00	1.00	1.00
  	1000	0	536	4.30	2.75	3.02	3.08
  		1		0.96	1.09	1.13	1.15
  		2		0.21	0.8	1.03	1.03
  		3		0.09	1.00	1.00	1.00
  		4		0.04	1.00	1.00	1.00
  0.7	40	0	997	3.27	2.15	2.09	2.11
  		1		0.85	1.02	1.09	1.10
  		2		0.21	0.98	1.02	1.02
  		3		0.06	0.99	0.99	0.99
  		4		0.01	1.00	1.00	1.00
  	500	0	1589	4.13	2.22	2.35	2.39
  		1		1.04	1.19	1.21	1.20
  		2		0.30	0.97	1.05	1.05
  		3		0.14	1.00	1.00	1.00
  		4		0.08	1.00	1.00	1.00
  	1000	0	1751	5.17	3.71	4.03	4.09
  		1		1.01	1.15	1.24	1.25
  		2		0.39	1.04	1.07	1.06
  		3		0.16	0.99	1.00	1.00
  		4		0.11	1.00	1.00	1.00

  . Similarly, both shrinkage estimators RPS and RSE outperforms the full ridge model estimator irrespective of the corrected sub-model selected. This is consistent with the asymptotic theory presented in Section 3.
• ${\Delta }^{*}$ which measures the degree of deviation from the Assumption 1 on the parameter space, it is clear that one cannot go wrong with the use of shrinkage estimators even if the selected sub-model is wrongly specified. As evident from Table 1 and Table 2, Figure 1 and Figure 2, if the selected sub-model is correct, that is, ${\Delta }^{*}=0$, then the shrinkage estimators are relatively efficient compared with the ridge full model estimator. On the other hand, if the sub-model is misspecified, the gain slowly diminishes. However, in terms of risk, the shrinkage estimators are at least as good as the full ridge model estimator. Therefore, the use of shrinkage estimators makes sense in application when a sub-model cannot be correctly specified.
• The RMSE of the ridge-type estimators are an increasing function of the amount of multicollinearity. This indicates that the ridge-type estimators perform better than the classical estimator in the presence of multicollinearity among predictor variables.

4.2. Comparison with LASSO-Type Estimators

We compare our listed estimators with the LASSO and adaptive LASSO estimators. A 10-fold cross-validation is used for selecting the optimal value of the penalty parameters that minimizes the mean square errors for the LASSO-type estimators. The results for $\rho =0.3,0.7,0.9$, $n=60,100$, $p_1=10$ and $p_2=50,500,1000,2000$ are presented in

Table 3. RMSEs of estimators with respect to ${\widehat{\mathit{\beta }}}^{\mathrm{RFM}}$ when $\Delta =0$ for $p_1=10$.

n	$\mathit{\rho }$	${\mathit{p}}_2$	CNI	RSM	RPT	RSE	RPS	LASSO	aLASSO
60	0.3	50	35.64	3.31	2.25	1.82	1.95	1.23	1.28
		500	452.76	4.13	3.71	2.61	3.01	1.47	1.52
		1000	1265.34	5.02	4.28	4.61	4.78	1.96	2.15
		2000	4567.56	7.13	5.10	6.18	6.39	2.70	3.06
	0.7	50	61.34	3.52	3.05	2.51	2.55	1.14	1.21
		500	743.17	4.49	3.65	3.41	3.50	1.36	1.58
		1000	2350.89	5.84	4.11	4.32	4.61	1.68	1.95
		2000	6908.39	8.10	5.31	6.24	6.29	1.84	2.02
	0.9	50	120.21	4.21	3.61	3.34	3.35	1.10	1.05
		500	950.98	4.82	3.3.8	3.72	3.73	1.21	1.16
		1000	5892.51	6.35	4.10	5.01	5.13	1.42	1.31
		2000	8352.73	8.51	4.63	5.24	5.38	1.61	1.35
100	0.3	50	31.21	2.91	2.54	2.12	2.23	1.32	1.36
		500	356.64	3.75	3.31	2.84	2.92	1.54	1.61
		1000	975.32	4.25	2.53	3.42	3.61	1.92	2.06
		2000	2764.84	5.61	4.25	4.91	5.08	2.31	2.46
	0.7	50	52.79	3.18	2.61	2.30	2.37	1.28	1.53
		500	578.43	4.28	3.05	3.52	3.59	1.46	2.07
		1000	1281.66	5.10	3.26	3.78	3.82	1.84	2.52
		2000	3498.30	6.12	3.01	4.26	4.33	2.27	2.41
	0.9	50	79.41	4.11	3.41	3.21	3.28	1.28	1.21
		500	681.43	4.35	3.55	3.41	3.50	1.43	1.51
		1000	1470.32	5.82	3.18	4.01	4.14	1.72	1.79
		2000	4105.90	7.04	4.57	5.22	5.32	1.87	1.96

. We observe the following from Table 3.

• The performance of the sub-model estimator is the best among all estimators.
• The pretest ridge estimator performs better than the other estimators. However, for larger values of sparse predictors $p_2$ the shrinkage estimators outperform the pretest estimator.
• The performance of the LASSO and aLASSO estimators are comparable when $\rho$ is small. The pretest and shrinkage estimators remain stable for a given value of $\rho$.
• For a large number of sparse predictors $p_2$, the shrinkage and pretest ridge estimators outperforms the lasso-type estimators. This indicates the superiority of the shrinkage estimators over the LASSO-type estimators. Therefore shrinkage estimators are preferable when there is multicollinearity in our predictor variables.

5. Real Data Application

We consider two real data analyses using Amsterdam Growth and Health Data and a genetic and brain network connectivity edge weight data to illustrate the performance of the proposed estimators.

5.1. Amsterdam Growth and Health Data (AGHD)

The AGHD data is obtained from the Amsterdam Growth and Health Study [25]. The goal of this study is to investigate the relationship between lifestyle and health in adolescence into young adulthood. The response variable Y is the total serum cholesterol measured over six time points. There are five covariates: $X_1$ is the baseline fitness level measured as the maximum oxygen uptake on a treadmill, $X_2$ is the amount of body fat estimated by the sum of the thickness of four skinfolds, $X_3$ is a smoking indicator (0 = no, 1 = yes), $X_4$ is the gender (1 = female, 2 = male), and time measurement as $X_5$ and subject specific random effects.

A total of 147 subjects participated in the study where all variables were measured at $n_i=6$ time occasions. In order to apply the proposed methods, firstly, we apply a variable selection based on AIC procedure to select the sub-model. For the AGHD data, we fit a linear mixed model with all the five covariates for both fixed and subject specific random effects by two stage selection procedure for the purpose of choosing both the random and fixed effects. The analysis found $X_2$ and $X_5$ to be significant covariates for prediction of the response variable serum cholestrol and the other variables are ignored since they are not significantly important. Based on this information, a sub-model is chosen to be $X_2$ and $X_5$ and the full model includes all the covariates. We construct the shrinkage estimators from the full-model and sub-model. In terms of null hypothesis, the restriction can be written as ${\mathit{\beta }}_2=({\beta }_1,{\beta }_3,{\beta }_4)=(0,0,0)$ with $p=5$, $p_1=2$ and $p_2=3$.

To evaluate the performance of the estimators, we obtain the mean square prediction error (MSPE) using bootstrap samples. We draw 1000 bootstrap samples of the 147 subjects from the data matrix $\{(Y_{ij},{\mathbf{X}}_{ij}),i=1,2,\dots ,147;j=1,2,\dots ,6\}$. We then calculate the relative prediction error (RPE) of ${\mathit{\beta }}_1^{*}$ with respect to ${\mathit{\beta }}_1^{\mathrm{RFM}}$, the full model estimator. The RPE is defined as

$$\mathrm{RPE}({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}:{\widehat{\mathit{\beta }}}_1^{*})={\displaystyle \frac{\mathrm{MSPE}\left({\widehat{\mathit{\beta }}}_1^{*}\right)}{\mathrm{MSPE}\left({\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}\right)}}={\displaystyle \frac{{(\mathbf{Y}-{\mathbf{X}}_1{\widehat{\mathit{\beta }}}_1^{*})}^{\prime }(\mathbf{Y}-{\mathbf{X}}_1{\widehat{\mathit{\beta }}}_1^{*})}{{(\mathbf{Y}-{\mathbf{X}}_1{\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}})}^{\prime }(\mathbf{Y}-{\mathbf{X}}_1{\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}})}},$$

where ${\mathit{\beta }}_1^{*}$ is one of the listed estimators. If $\mathrm{RPE}<1$, then ${\widehat{\mathit{\beta }}}_1^{*}$ outperforms ${\widehat{\mathit{\beta }}}_1^{\mathrm{RFM}}.$

Table 4. Estimate, standard error for the active predictors and RPEs of estimators with respect to full-model estimator for the Amsterdam Growth and Health Study data.

	RFM	RSM	RPT	RSE	RPS	LASSO	aLASSO
Estimate(${\mathit{\beta }}_2$ )	0.381	0.395	0.392	0.389	0.390	0.624	0.611
Standard error	0.104	0.102	0.100	0.009	0.008	0.081	0.079
Estimate (${\mathit{\beta }}_5$)	0.137	0.125	0.131	0.130	0.133	0.101	0.105
Standard error	0.012	0.010	0.009	0.011	0.010	0.013	0.012
RPE	1.000	0.723	0.841	0.838	0.831	0.986	0.973

reports the estimates, standard error of the non-sparse predictors and RPEs of the estimators with respect to the full model. As expected, the sub-model ridge estimator ${\widehat{\mathit{\beta }}}_1^{\mathrm{RSM}}$ has the minimum RPE because it is computed when the sub-model is correct, that is, ${\Delta }^{*}=0$. It is evident by the RPE values in Table 4 that the shrinkage estimators are superior to the LASSO-type estimators. Furthermore, the positive shrinkage is more efficient than the shrinkage ridge estimator.

5.2. Resting-State Effective Brain Connectivity and Genetic Data

This data comprises longitudinal resting-state functional magnetic resonance imaging (rs-fMRI) effective brain connectivity network and genetic study [26] data obtained from a sample of 111 subjects with a total of 319 rs-fMRI scans from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database. The 111 subjects comprise 36 cognitively normal (CN), 63 mild cognitive impairment (MCI) and 12 Alzheimer’s Disease (AD) subjects. The response is a network connection between regions of interest estimated from an rs-fMRI scan within the Default Mode Network (DMN), and we observe a longitudinal sequence of such connections for each subject with the number of repeated measurements. The DMN consists of a set of brain regions that tend to be active in resting-state, when a subject is mind wandering with no intended task. For this data analysis, we consider the network edge weight from the left intraparietal cortex to posterior cingulate cortex (LIPC → PCC) as our response. The genetic data are single nucleotide polymorphism (SNPs) from non-sex chromosomes, i.e., chromosome 1 to chromosome 22. SNPs with minor allele frequency less than 5% are removed as are SNPs with a Hardy–Weinberg equilibrium p-value lower than ${10}^{-6}$ or a missing rate greater than 5%. After preprocessing we are left with 1,220,955 SNPs and the longitudinal rs-fMRI effective connectivity network using the 111 subjects with rs-fMRI data. The response is network edge weight. There are SNPs which are the fixed effects and subject specific random effects.

In order to apply the proposed methods, we use a genome- wide association study (GWAS) for screening the genetic data to 100 SNPs. We implement a second screening by applying multinomial logistic regression to identify a smaller subset of the 100 SNPs that are potentially associated with disease (CN/MCI/AD). This yields a subset of top 10 SNPs. This showed the top 10 SNPs are the most important predictors and the other 90 SNPs are ignored as not significant. We now have two models, which are the full model with all 100 SNPs and sub-model with 10 SNPs selected. Finally, we construct the pretest and shrinkage estimators from the full-model and sub-model.

We draw 1000 bootstrap samples with replacements from the corresponding data matrix $\{(Y_{ij},{\mathbf{X}}_{ij}),i=1,\dots ,111;j=1,\dots ,n_i\}.$ We report the RPE of the estimators based on the bootstrap simulation with respect to the full model ridge estimator in

Table 5. RPEs of estimators.

	RFM	RSM	RPT	RSE	RPS	LASSO	aLASSO
RPE	1.000	0.802	0.947	0.932	0.928	1.051	1.190

. We observe that the RPE of the sub-model, pretest, shrinkage and positive shrinkage ridge estimators outperforms the full model estimator. Clearly, the sub-model ridge estimator has the smallest RPE since it’s computed when the candidate sub-model is correct, i.e., $\Delta =0$. Both shrinkage ridge estimators outperform the pretest ridge estimator. Particularly, the positive shrinkage performed better than the shrinkage estimator. The performance of both shrinkage and pretest ridge estimators are better than the LASSO-type estimators. Thus, the data analysis is in line with our simulation and theoretical findings.

6. Conclusions

In this paper, we present efficient estimation strategies for the linear mixed effect model when there exists multicollinearity among predictor variables for high-dimensional data application. We considered the estimation of fixed effects parameters in the linear mixed model when some of the predictors may have a very weak influence on the response of interest. We introduced pretest and shrinkage estimation in our model using the ridge estimation as the reference estimator. In addition, we established the asymptotic properties of the pretest and shrinkage ridge estimators. Our theoretical findings demonstrate that the shrinkage ridge estimators outperform the full model ridge estimator and perform relatively better than the sub-model estimator in a wide range of the parameter space.

Additionally, a Monte Carlo simulation was conducted to investigate and assess the finite sample behavior of proposed estimators when the model is sparse (restrictions on parameters hold). As expected, the sub-model ridge estimator outshines all other estimators when the restrictions hold. However, when this assumption is violated, the shrinkage and pretest ridge estimators outperform the sub-model estimator. Furthermore, when the number of sparse predictors are extremely large relative to the sample size, the shrinkage estimators outperform the pretest ridge estimator. These numerical results are consistent with our asymptotic result. We also assess the relative performance of the LASSO-type estimators with our ridge-type estimators. We observe that the performance of pretest and shrinkage ridge estimators are superior to the LASSO-type estimators when predictors are highly correlated. For our real data application, the shrinkage ridge estimators are superior with the smallest relative prediction error compared to the LASSO-type estimators.

In summary, the results of the data analyses strongly confirm the findings of the simulation study and suggest the use of the shrinkage ridge estimation strategy when no prior information about the parameter subspace is available. The results of our simulation study and real data application are consistent with available results in [27,28,29].

In our future work, we will focus on other penalty estimators like the Elastic-Net, the minimax concave penalty (MCP), and the smoothly clipped absolute deviation method (SCAD) as estimation strategy in LMM for high-dimensional data. These estimators will be assessed and compared with the proposed ridge-type estimators. Another interesting extension will be integrating two sub-models by incorporating ridge-type estimation strategies in the linear mixed effect models. The goal is to improve the estimation accuracy of the non-sparse set of the fixed effects parameters by combining an over-fitted model estimator with an under-fitted one [27,29]. This approach will include combining two sub-models produced by two different variable selection techniques from the LMM [28].