Small area quantile estimation via spline regression and empirical likelihood

Section 5. Monte Carlo simulations

• Previous
• Next

In this section, we use simulation to evaluate the performances of the proposed penalized spline regression model based empirical likelihood estimators (PEL) and their MSE estimates. When only the covariate population means are known the proposed estimators are compared with only the nested-error linear regression model based empirical likelihood estimator (LEL) of Chen and Liu (2018), and the direct estimator (DE). When covariate values are known for all sample units, the comparison is extended to also include six estimators of Tzavidis et al. (2010), denoted as EBLUP/naïve, EBLUP/CD, EBLUP/RKM, M-quantile/naïve, M-quantile/CD and M-quantile/RKM. Here, EBLUP/CD and M-quantile/CD denote the EBLUP and M-quantile estimator are obtained based on the CDF proposed by Chambers and Dunstan (1986), and corresponding estimators based on the CDF proposed by Rao, Kovar and Mantel (1990) write as RKM.

Similar to Chen and Liu (2018), we must choose $q\left(u\right)$ in the DRM. Two candidates are $q_1\left(u\right)\mathrm{<mtext> </mtext>\; <mo>=</mo>\; <mtext> </mtext>}{\left(\mathrm{1,}u\right)}^{\prime }$ and $q_2\left(u\right)\mathrm{<mtext> </mtext>\; <mo>=</mo>\; <mtext> </mtext>}{\left(\mathrm{1,}\text{sign}\left(u\right)\sqrt{\left|u\right|}\right)}^{\prime }.$ Some preliminary simulation results indicate that $q_1\left(u\right)\mathrm{<mo>=</mo>}{\left(\mathrm{1,}u\right)}^{\prime }$ works well for the P-splines fitted non-parametric regression model, but $q_2\left(u\right)$ does not. Instead, the choice of $q_2^{*}\left(u\right)\mathrm{<mo>=</mo>}{\left(\mathrm{1,}u\mathrm{,}{u}^2\right)}^{\prime }$ leads to competitive performance. So, we use $q_1\left(u\right)$ and $q_2^{*}\left(u\right)$ in our simulation.

Following Rao et al. (2014) and Torabi and Shokoohi (2015), we generated data from the following three models:

$$\begin{array}{ll}\text{A}:\hfill & y_{ij}\mathrm{<mo>=</mo>\; 1}+x_{ij}+v_i+{\epsilon }_{ij}\mathrm{,}\hfill \\ \text{B}:\hfill & y_{ij}\mathrm{<mo>=</mo>\; 1}+x_{ij}+x_{ij}^2+v_i+{\epsilon }_{ij}\mathrm{,}\hfill \\ \text{C}:\hfill & y_{ij}\mathrm{<mo>=</mo>\; 1}-x_{ij}+\text{0}\text{.5}\exp \left(x_{ij}\right)+v_i+{\epsilon }_{ij}\mathrm{.}\hfill \end{array}$$

They lead to linear, quadratic and exponential regression functions respectively. We set the number of small areas to be 30 and area population sizes $N_i\mathrm{<mo>=</mo>\; 500}\left(i+1\right)\mathrm{,}i\mathrm{<mo>=</mo>\; 0,}\mathrm{1,}\dots \mathrm{,}29.$ We generated covariate $x_{ij}$ from $N\left(\mathrm{0,}1\right).$ Once $x_{ij}$ are generated, we treated them as fixed in the simulation. The area-specific random effect $v_i$ were generated from $N\left(\mathrm{0,}1\right),$ and the errors ${\epsilon }_{ij}$ were generated from the following four distributions.

$$\begin{array}{ll}\left(\text{i}\right)\hfill & :N\left(\mathrm{0,}1\right)\mathrm{,}\hfill \\ \left(\text{ii}\right)\hfill & :t\left(3\right)\mathrm{,}\hfill \\ \left(\text{iii}\right)\hfill & \text{:}\text{normal}\text{mixture }\text{0}\text{.5}N\left(-\mathrm{1,}1\right)+\text{0}\text{.5}N\left(\mathrm{1,}1\right)\mathrm{,}\hfill \\ \left(\text{iv}\right)\hfill & :N\left(\mathrm{0,}{\sigma }_i^2\right)\mathrm{,}\text{with }{\sigma }_i\sim U\left(\text{0}\text{.5}\mathrm{,}2\right)\mathrm{,}i\mathrm{<mo>=</mo>\; 0,}\dots \mathrm{,}29.\hfill \end{array}$$

Distribution (ii) has a heavy tail, distributions (ii) and (iii) are symmetric, and distribution (iv) is heteroscedastic.

We used $R\mathrm{<mo>=</mo>}\text{1,000}$ repetitions in the simulation and drew random samples of size $n\mathrm{<mo>=</mo>\; 500}$ without replacement from the population in each repetition. To avoid the possibility that some small areas have too few sample units, we drew $n-60$ units at the population level and allocated an additional 2 units in each small area. We used R package mgcv for the REML method with default options for values of $p$ and $K$ when fitting the P-spline function (2.4). We calculated estimates of the 5%, 25%, 50%, 75%, and 95% small area quantiles denoted as DE, LEL1, LEL2, PEL1, PEL2, for direct estimator, estimators of Chen and Liu (2018) and the proposed estimators using $q_1(\cdot )$ and $q_2(\cdot ).$ We report their average mean squared error (AMSE) and absolute biases (ABIAS) defined below:

$$\begin{array}{ll}\text{AMSE}\hfill & \mathrm{<mo>=</mo><mtext> </mtext>}{\left\{R\left(m+1\right)\right\}}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 0}}^m}{\displaystyle \sum _{r\mathrm{<mo>=</mo>\; 1}}^R{\left({\widehat{\xi }}_i^{\left(r\right)}-{\xi }_i^{\left(r\right)}\right)}^2}\mathrm{,}\hfill \\ \text{ABIAS}\hfill & \mathrm{<mo>=</mo>}{\left(m+1\right)}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 0}}^m}\left|R^{-1}{\displaystyle \sum _{r\mathrm{<mo>=</mo>\; 1}}^R{\widehat{\xi }}_i^{\left(r\right)}}-R^{-1}{\displaystyle \sum _{r\mathrm{<mo>=</mo>\; 1}}^R{\xi }_i^{\left(r\right)}}\right|\mathrm{,}\hfill \end{array}$$

where ${\widehat{\xi }}_i^{\left(r\right)}$ is either one of the quantile estimates of for the $i^{\text{th}}$ small area in the $r^{\text{th}}$ repetition. The results under Models A, B, and C are given in Tables 5.1-5.3 respectively. Both PEL and LEL are based on ${\widehat{F}}_i^{\left(a\right)}$ and its mirror version in Chen and Liu (2018).

Under Model A, the linear model is valid. Hence, we expect LEL to be superior. According to Table 5.1, two methods are similar for the 25%, 50% and 75% quantiles. LELs outperform PELs for the 5% quantile while the comparison reverses for the 95% quantile. Both PEL and LEL outperform DE for the 25%, 50% and 75% quantiles with big margins. An overall impression is that the proposed methods still work satisfactorily.

Under Model B, the linear model breaks down mildly. Results in Table 5.2 show that the PEL estimators have lower AMSE for lower quantiles. The LELs still have low AMSE in spite of have higher ABIAS. The advantage of the proposed PEL under the non-parametric nested-error regression models focus for quantiles in middle levels. With fewer observations near extreme quantiles, the non-parametric model is hard to fit.

The linearity is seriously violated under Model C. LEL is expected to have poor performance and this is evident as shown in Table 5.3. At the same time, PELs work well for the 25%, 50% and 75% quantiles. The choice of $q_2^{*}\left(u\right)$ also helps in general. For extreme quantiles, PELs remain unworth the trouble compared with DE.

Next, we study estimators applicable when covariate values are known for all sample units. The simulation includes EB0, EB1, EB2, MQ0, MQ1 and MQ2 stand for EBLUP/naïve, EBLUP/CD, EBLUP/RKM, M-quantile/naïve, M-quantile/CD and M-quantile/RKM respectively. We set relatively small population sizes $N_i\mathrm{<mo>=</mo>\; 500}$ to save some computation. Table 5.4 contains the AMSE of these estimators under Models A, B and C with $N\left(\mathrm{0,}1\right)$ error distribution. To save space, we do not present the corresponding bias results. The simulation results show that the proposed method has lower AMSE and ABIAS (not presented) in general. It works well even for quantiles at rather extreme levels.

To save space, we pool the AMSE results for all 5 levels of quantiles in Table 5.5. The entry corresponding to $A_i$ is the average AMSE for estimating quantiles at levels 5%, 25%, 50%, 75%, and 95% when data are generated from Model A with error distribution (i). We notice that with more detailed information on covariates, the LEL and PEL estimators are substantially more accurate compared to results in Tables 5.1-5.3. From Model A to Model C, the regression line becomes less linear. Correspondingly, the proposed quantile estimators have greater advantages against other estimators.

Now we evaluate the bootstrap MSE estimator proposed in Section 4. Because this method involves heavy computation, we confined the simulation to the estimator based on ${\widehat{F}}_i^{\left(b\right)}\left(u\right)$ with basis function $q_1\left(u\right)\mathrm{<mo>=</mo>}{\left(\mathrm{1,}u\right)}^{\prime }$ and put $B\mathrm{<mo>=</mo>\; 100,}L\mathrm{<mo>=</mo>\; 100.}$ We report the average ratios of the estimated MSEs and the simulated MSEs across all the small areas. The closer the ratio to one, more accurate the bootstrap MSE estimate. From Table 5.6 we can see that the average ratios close to one in majority situations except for error distribution (iv) on extreme levels of quantiles. We conclude that the bootstrap MSE estimator is generally satisfactory.

• Previous
• Article main page
• Next