A comparative study of estimators in multilevel linear models

Multilevel Models are widely used in organizational research, educational research, epidemiology, psychology, biology and medical fields. In this paper, we recommend the situations where Bootstrap procedures through Minimum Norm Quadratic Unbiased Estimator (MINQUE) can be extremely handy than that of Restricted Maximum Likelihood (REML) in multilevel level linear regression models. In our simulation study the bootstrap by means of MINQUE is superior to REML in conditions where normality does not hold. Moreover, the real data application also supports our findings in terms of accuracy of estimates and their standard errors.


Introduction
Multilevel data or clustered data are commonly observed in schools, health institutions, and epidemiology. Multilevel models are also called hierarchical, mixed effects, or random effects models Snijders and Bosker [1], Raudenbush and Bryk [2].
Maximum likelihood (ML) method estimates and estimates standard errors were used by Maas and Hox [3]. Wen et al. [4] concluded that Bayesian spatial-temporal model is superior to the random effects model and spatial model for investigating the effects of weather and roadway characteristics on crash incidence. Brown and Draper [5] utilized ML method of estimation and accomplished that in small sample sizes the estimates are biased. MINQUE recommended by Rao [6], as an alternate to ML estimator. The method, however, does not rely on the assumption of normality in multilevel linear models. According to Bagakas [7], one major problem with the MINQUE estimators is that standard errors of the minimum norm quadratic unbiased estimators cannot be computed because of the non-existence of formulae. In situations, where a researcher attempts to construct confidence interval and perform testing of hypothesis about the parameter then the MINQUE is not appropriate. The researcher then needs to use an alternate scheme such as bootstrapping, where not only the parameter estimates but also their standard errors can be estimated by applying different estimation methods such as MINQUE or ML method of estimation.
In practice, both parametric and nonparametric bootstrap can be used. However, when the assumption of normality does not exist the nonparametric bootstrap is handy. As the MINQUE method of estimation is free from the normality assumption, so the bootstrap by means of MIN-QUE will be used. Swallow  Bagakas [7] used bootstrap by means of MINQUE. Similarly, Meijer et al. [9] concluded that multilevel bootstrapping performance was excellent in small sample sizes in multilevel models. Carpenter et al. [10] carried out a simulation study where they compared the relative performance of parametric bootstrap and nonparametric residuals bootstrap methods by using multilevel linear models. Hutchison et al. [11] successfully carried out simulation study on a two-level model. They applied the procedure of nonparametric cases bootstrap and promising standard errors of the estimates were obtained. Wang et al. [12] used multilevel linear model to apply nonparametric residual bootstrap through a SAS macro. Nonparametric residual bootstrap estimates standard errors were promising. Delpish [13] also compared REML and Bootstrap by means of MINQUE in her study. Ali et al., [14] concluded that ML gave better results than Penalized Quasilikelihood (PQL) for small sample conditions in multilevel model. To get accurate estimates of both fixed and random effects ML requires relatively small sample compared to PQL in multilevel logistic models (Ali et al. [15]). In a study by Zeng et al. [16] revealed that univariate spatial model gave lower deviance information criteria (DIC) and accurate estimates of parameters as compared to bivariate spatial model while investigating the factors responsible for vehicle crash on freeway. The proposed multivariate random-parameters spatio-temporal Tobit model gave lower Deviance Information Criteria (DIC), Mean Absolute Deviance (MAD) and Mean Squared Prediction Errors (MSPE) then the competing model such as multivariate random-parameters Tobit model and a multivariate random-parameters spatial Tobit model (Zeng et al. [17]. It was confirmed from the results that spatio-temporal correlation and interaction have significance in the area wide crash data.
In this paper, the researchers investigate the performance of REML and Bootstrap by means of MINQUE under varying conditions of the number of groups, Intra-class correlation and different skewed distributions.

Materials and methods
For this study a random intercept and random slope multilevel linear model was used. The model has single explanatory variable at each level. The model is given below: The combined model was obtained by substituting level 2 model in level 1 model: (Fixed part)+(Random part) Where X ij is the Level 1 explanatory variable, W j corresponds to Level 2 explanatory variable, γ 00 , γ 10 , γ 01 and γ 11 are the fixed effects, e ij is assumed to follows a normal distribution i.e e ij * N (0, s 2 e ). In case of normality, u oj and u 1j assumed to follow a multivariate normal distribution as PLOS ONE s 2 u Corresponds to the random intercept variance, s 2 1 is the random slope variance and σ u1 is the covariance term.

Design factors
1. Three levels of number of groups were used in this study: 30,100 and 120.
2. Three levels of intra-class correlations were used: 0.01, 0.10 and 0.20. Where the intra-class correlation coefficient (ICC) is given as 3. Three distributions were used: Normal distribution, Lognormal distribution and Exponential distribution

Analysis
Two estimation procedures Restricted Maximum Likelihood and Bootstrap by means of MIN-QUE were used in all the three distribution conditions. All the simulations and bootstrapping were performed in SAS 9.2 to obtain estimates and their standard errors. Algorithm. The procedure for cases bootstrap is as given below: 1. Draw with replacement J group level units along with corresponding scores on group level variable W � j .
2. Then draw with replacement n j individual level units within group level unit j, j = 1, 2. . .. . .. . ., J. This results the bootstrap data (Y � , X � ) and this data set is then combined with the group level variable W � j in order to get (Y � , X � , W � j ) the desired bootstrap sample.
3. Obtain the minimum norm quadratic unbiased estimates of the model parameters from the bootstrap replicated sample.

Obtain the mean value of estimates by usinĝ
And the bootstrap parameter estimate standard error is obtain as The real data was selected from High School & Beyond Survey data set, which is a national survey of United States conducted by National Center for Educations Statistics (NCES) about Public and Catholic schools. For the purpose of illustration, a dataset of 30 schools was randomly selected from the data of 160 schools.

Results
Tables 1 and 2 show that the bootstrap procedure showed perfect results in terms of accuracy of the fixed and random effects estimates, however, REML method estimates were comparable to that of the bootstrap procedure at 100 and 1200 groups respectively. Similarly, from Table 3 it is evident that the bootstrap CI outclassed the REML CI at the first two levels of the number of group (30 and 100) factor when the distribution was normal.
The bootstrap procedure was superior to REML in terms of accuracy of the fixed effects and random effect estimates as can be seen in Tables 4 and 5 for lognormal distribution. Moreover, Table 6 reveals that the bootstrap CI outperformed the REML CI at all levels of the number of groups when data was generated from lognormal distribution. Furthermore, when the distribution of the data was exponential again the bootstrap method outshined the REML method estimates as shown in Tables 7-9 respectively.

Real data application
The application of bootstrap by means of MINQUE method to the real data is demonstrated in this section. A two-level model was fitted to a subsample data drawn from High School & Beyond (HSB) data. The data consist of two levels i.e school level and student level. HSB data consists of 7185 students nested within 160 schools. The data contains four level 1 or individual level variables and six level 2 or group level variables in total. For the purpose of illustration of bootstrap by means of MINQUE method only 30 schools were drawn randomly from 160 schools. The total numbers of level 1 units are 1447 and level 2 units are 30. Students MATH Table 1. Average relative parameter bias of fixed effects estimates obtained for normal distribution data (First = REML estimation procedure, Second = Bootstrap estimates are enclosed in parenthesis). ACHIEVEMENT SCORE was taken as a response variable, SES was selected as a level 1 variable and MEANSES was selected as a level 2 variable. A two-level model used in this real data application is given below

Groups
Level-1 model Level-2 models . . . n j and j ¼ 1; 2; . . . J The combined model can be written as REML and bootstrap by means of MINQUE estimation procedures were used to estimate both fixed effects and random effects using HSB: 30 schools data set for the model in equation (1.8). The SAS package procedure PROC MIXED was used to obtain REML estimates and estimates standard errors. The REML confidence intervals were then constructed for each   8), B = 1000 bootstrap replicates were obtained using cases bootstrap. The mean of 1000 estimates were then taken to obtain the bootstrap estimate. This means that the bootstrap estimate of any parameter is the average of one thousand estimates. On the other hand, single estimate for each parameter was obtained under REML method of estimation. Bootstrap confidence intervals were constructed for each parameter in the model by using the percentile method. The data set of 30 schools randomly selected from 160 school's data is presented in Table 10. Table 11 illustrates estimates and estimated standard errors under REML and bootstrap by means of MINQUE methods of estimation. Moreover, 95% CI's are also given in Table 11. There is not much difference to choose between the two procedures as for as the accuracy of the estimates is concerned. However, both fixed and random effects estimate standard errors were lower under bootstrap by means of MINQUE. The widths of the REML CI's were clearly higher than that of the percentile bootstrap CI's. Overall, for real data, bootstrap by means of MINQUE performs better than that of the REML method of estimation especially in terms of precision. Simulation results also exposed that bootstrap by means of MINQUE procedure outperformed the REML method of estimation particularly in terms of estimates promising standard errors.

Conclusion
REML produced unbiased fixed effects estimates at the second level and third level of the number of groups (100 and 120) factor. On the other hand, the bootstrap fixed effects estimates Table 5. Average relative parameter bias of the random effect estimates obtained for lognormal distribution data (First = REML estimation procedure, Second = Bootstrap estimates are enclosed in parenthesis).  were unbiased across all conditions. Additionally, the bootstrap procedure outperformed the REML method in terms of accuracy of the random effects estimates when the number of groups was 30. Based on the above normal data results, it is recommended that at least 30 groups are essential to obtain unbiased fixed effects estimates and their standard errors under REML method of estimation. Furthermore, 100 groups are essential to achieve accurate random effects estimates and their standard errors under REML method of estimation. It is also   recommended that bootstrap by means of MINQUE can be superior to REML when the number of groups are 30 and normality holds.

Groups
In general, the estimates and estimated standard errors were biased for the two skewed distribution data when the number of groups was 30 under REML method of estimation. On the other hand, the bootstrap estimates and estimated standard errors were unbiased across all conditions. To put it differently, the bootstrap fixed effects and random effects estimates coverage rates were not only acceptable but also superior to that of REML estimates coverage rates across all conditions. Furthermore, REML level 2 random effects estimates coverage rates were It is recommended on the basis of these study results, whenever the data are based on skewed distributions and normality assumption does not hold, REML should not be used particularly for inference. In such situations, the bootstrap standard errors by means of MINQUE can be used for inference to achieve precise results.