6. Simulation studies

Cyril Favre Martinoz, David Haziza and Jean-François Beaumont

Previous | Next

6.1  Winsorization in a simple random sampling without-replacement design

We carried out a simulation study to examine the properties of several robust estimators using 11 populations. The first 10 of size $N=\text{5,000}$ consists of a variable of interest $y.$ In each population, the $y-$ values were generated according to the following model:

$$Y_i=U_i+{\delta }_i{V}_i\mathrm{,}$$

where $U_i,{\delta }_i$ and $V_i$ are random variables whose distributions are described in Table 6.1. Population 1 was generated according to a normal distribution. Populations 2 through 5 were generated using a mixture of normal distributions with contamination rates ranging from 0.5% to 5%. Populations 6 through 8 were generated according to skewed distributions. Populations 9 and 10 were generated using a mixture of lognormal distributions with contamination rates equal to 0.5% and 5%. Population 11 of size $N=\text{5,000}$ is from the information technology survey produced by the French National Institute for Statistics and Economic Studies (INSEE) in 2011. One of the survey's objectives is to estimate the e-commerce sales of French companies. We use the "sales� variable in our simulation. The distribution of $y$ in each population is plotted in Figure 6.1. In addition, Table 6.2 presents a number of descriptive statistics for each of the populations used. For confidentiality reasons, the units for Population 11 are not shown in the plot. Similarly, there are no descriptive statistics for Population 11 in Table 6.2.

In each population, we selected $M=\text{5,000}$ samples according to a simple random sampling without-replacement design of size $n=100,300$ and $500.$ For each sample, we calculated the expansion estimator $\widehat{t}$ and the robust estimator (4.8). Let $y_{\left(1\right)}\mathrm{,}\dots \mathrm{,}{y}_{\left(n\right)}$ be the values of the $y-$ variable arranged in ascending order. We also calculated the first-, second- and third-order winsorized estimators, where the $p^{\text{th}}-$ order winsorized estimator is obtained by replacing the $p$ largest values in the sample with the value $y_{\left(n-p\right)}\mathrm{,}p=1,2,3.$ In a classical statistical context, Rivest (1994) showed that the first-order winsorized estimator has good mean-square-error properties for a large class of skewed distributions.

As a measure of the bias of an estimator $\widehat{\theta },$ we calculated the Monte Carlo relative bias (in percentage):

$${\text{BR}}_{\text{MC}}\left(\widehat{\theta }\right)=\frac{\frac{1}{M}{\displaystyle \sum _{m=1}^M}\left({\widehat{\theta }}_{(m)}-t\right)}{t}\times \mathrm{100,}$$

where ${\widehat{\theta }}_{\left(m\right)}$ denotes the estimator $\widehat{\theta }$ in sample $m\mathrm{,}m=1,\dots \mathrm{,}5,000.$ We also calculated the relative efficiency of the robust estimators with respect to the expansion estimator, $\widehat{t}:$

$${\text{RE}}_{\text{MC}}\left(\widehat{\theta }\right)=\frac{\frac{1}{M}{\displaystyle \sum _{m=1}^M}{\left({\widehat{\theta }}_{(m)}-t\right)}^2}{\frac{1}{M}{\displaystyle \sum _{m=1}^M}{\left({\widehat{t}}_{(m)}-t\right)}^2}\times 100.$$

The results are shown in Table 6.3.

The results presented in Table 6.3 show that the once-winsorized estimator has lower bias and is generally more efficient than the two times and three times winsorized estimators, which is consistent with the results obtained by Rivest (1994). It is interesting to compare the robust estimator ${\widehat{t}}_R$ and the once-winsorized estimator. In the case of Population 1, which does not contain any influential values, we see that both estimators have low bias and are as efficient as the expansion estimator. In the case of the populations with a mixture of normal distributions (Populations 2 to 5), we observe that the once-winsorized estimator is less efficient than the robust estimator in every scenario except for Population 5 with $n=300.$ In fact, the once-winsorized estimator is less efficient than the expansion estimator in every scenario except for Population 2 with $n=100.$ The robust estimator is more efficient than the expansion estimator except in Populations 4 and 5, for which we observe values of relative efficiency ranging from 91% to 102%. In the case of the populations with a mixture of lognormal distributions (Populations 9 and 10), we see that the bias and efficiency performance of the once-winsorized estimator and the robust estimator is very similar in all scenarios. The same is true for the skewed populations (Populations 6 to 8), for which the two estimators produce similar results. In the case of Population 11, the robust estimator has a lower bias than the once-winsorized estimator for $n=100,$ though it is less efficient (41% versus 47%). For $n=300$ and $n=500,$ the robust estimator has a lower bias and is significantly more efficient than the once-winsorized estimator.

Table 6.1
Models used to generate the populations
Table summary
This table displays the results of Models used to generate the populations. The information is grouped by Population (appearing as row headers), $U_i$  distribution, Mixture, ${\delta }_i$  distribution and $V_i$   distribution (appearing as column headers).

Population	$U_i$  distribution	Mixture	${\delta }_i$  distribution	$V_i$  distribution
1	$\mathcal{N}\left(\text{2,000; 500}\right)$	No
2	$\mathcal{N}\left(\text{2,000; 500}\right)$	Yes	$ℬ\left(\text{0}\text{.005}\right)$	$\mathcal{N}\left(\text{50,000; 10,000}\right)$
3	$\mathcal{N}\left(\text{2,000; 500}\right)$	Yes	$ℬ\left(\text{0}\text{.01}\right)$	$\mathcal{N}\left(\text{50,000; 10,000}\right)$
4	$\mathcal{N}\left(\text{2,000; 500}\right)$	Yes	$ℬ\left(\text{0}\text{.02}\right)$	$\mathcal{N}\left(\text{50,000; 10,000}\right)$
5	$\mathcal{N}\left(\text{2,000; 500}\right)$	Yes	$ℬ\left(\text{0}\text{.05}\right)$	$\mathcal{N}\left(\text{50,000; 10,000}\right)$
6	$ℒ\text{og}-\mathcal{N}\left(\text{log}\left(\text{2,000}\right);\text{1}\text{.2}\right)$	No
7	$ℒ\text{og}-\mathcal{N}\left(\text{log}\left(\text{2,000}\right);\text{1}\text{.5}\right)$	No
8	$ℱ\text{rechet}\left(\text{2,000};\text{2}\text{.5};\text{2}\text{.1}\right)$	No
9	$ℒ\text{og}-\mathcal{N}\left(\text{log}\left(\text{2,000}\right);\text{1}\text{.2}\right)$	Yes	$ℬ\left(\text{0}\text{.05}\right)$	$ℒ\text{og}-\mathcal{N}(\text{log}\left(\text{5,000}\right)\mathrm{;}\text{1}\text{.2)}$
10	$ℒ\text{og}-\mathcal{N}\left(\text{log}\left(\text{2,000}\right);\text{1}\text{.2}\right)$	Yes	$ℬ\left(\text{0}\text{.05}\right)$	$ℒ\text{og}-\mathcal{N}\left(\text{log}\left(\text{50,000}\right)\mathrm{;}\text{1}\text{.2}\right)$

 

Table 6.2
Descriptive statistics for the ten simulated populations
Table summary
This table displays the results of Descriptive statistics for the ten simulated populations. The information is grouped by Descriptive statistic (appearing as row headers), Population (appearing as column headers).

Descriptive statistic	Population
	1	2	3	4	5	6	7	8	9	10
min	132.3	314.9	105.3	275.9	187.4	23.6	7.6	2,000.9	20.5	26.6
max	3,968	79,506	78,526	80,540	78,690	252,612	379,751	2,159	305,612	1.3×106
$Q1$	1,639	1,667	1,664	1,666	1,685	883	743	200	920	913
Median	1,986	1,993	1,997	2,015	2,053	1,996	1,981	2,002	2,167	2,041
$Q3$	2,330	2,337	2,339	2,349	2,421	4,505	5,337	2,004	5,018	4,927
Mean	1,985	2,267	2,536	2,976	4,661	4,005	6,118	2,004	4,738	7,883
Standard deviation	503	3,709	5,506	7,119	11,470	7,353	17,190	5,89	9,796	33,111
Skewness	0.0	14.0	10.2	7.3	4.3	4.2	11.6	11.8	12.1	18.4
Kurtosis	3	209	109	56	20	19	196	228	267	570
CV	0.25	1.6	2.2	2.4	2.5	1.8	2.8	2.9×10-3	2.0	4.2

 

Figure 6.1 Distribution of the variable of interest in the 11 populations

Description for Figure 6.1

Table 6.3
Monte Carlo relative bias (in %) and relative efficiency (in parentheses) of several estimators
Table summary
This table displays the results of Monte Carlo relative bias (in %) and relative efficiency (in parentheses) of several estimators. The information is grouped by Population (appearing as row headers), $n$, ${\widehat{t}}_R$ and Winsorization (appearing as column headers).

Population	$n$	${\widehat{t}}_R$	Winsorization
			Once	Two times	Three times
1	100	-0.1(100)	-0.1(100)	-0.2(101)	-0.3(102)
	300	0.0(100)	-0.0(100)	-0.0(100)	-0.1(100)
	500	0.0(100)	-0.0(100)	-0.0(100)	-0.0(100)
2	100	-4.9(59)	-7.5(87)	-10.7(65)	-11.9(55)
	300	-2.9(87)	-3.0(129)	-6.8(158)	-9.5(169)
	500	-1.9(96)	-1.2(122)	-3.6(175)	-6.5(226)
3	100	-6.9(74)	-8.9(122)	-16.5(119)	-20.0(107)
	300	-3.5(99)	-1.9(122)	-5.6(171)	-10.6(232)
	500	-2.4(102)	-0.9(107)	-2.2(130)	-4.5(186)
4	100	-7.6(91)	-6.2(131)	-15.5(169)	-24.4(194)
	300	-2.9(101)	-0.6(103)	-2.1(118)	-4.4(154)
	500	-2.0(102)	-0.6(102)	-1.1(101)	-1.8(108)
5	100	-5.7(102)	-1.1(104)	-4.1(126)	-9.7(173)
	300	-2.2(102)	-0.4(100)	-0.8(101)	-1.4(102)
	500	-1.2(100)	-0.1(100)	-0.3(100)	-0.5(101)
6	100	-5.7(79)	-5.4(75)	-8.2(80)	-10.6(89)
	300	-2.6(84)	-2.6(79)	-3.9(81)	-5.1(88)
	500	-2.0(86)	-2.0(81)	-3.0(82)	-3.8(88)
7	100	-8.4(72)	-9.3(73)	-14.7(72)	-18.7(79)
	300	-4.5(86)	-4.4(95)	-7.8(91)	-10.2(95)
	500	-3.5(94)	-3.1(105)	-6.0(106)	-8.1(109)
8	100	-0.0(69)	-0.0(75)	-0.0(77)	-0.0(85)
	300	-0.0(82)	-0.0(88)	-0.0(87)	-0.0(95)
	500	-0.0(88)	-0.0(96)	-0.0(94)	-0.0(100)
9	100	-5.7(73)	-5.8(71)	-9.5(72)	-12.4(80)
	300	-3.5(87)	-3.5(85)	-5.4(88)	-6.8(98)
	500	-2.4(88)	-2.4(88)	-3.8(90)	-4.9(97)
10	100	-13.5(68)	-15.0(70)	-24.6(76)	-31.7(89)
	300	-7.5(80)	-7.2(79)	-12.1(85)	-16.3(97)
	500	-5.3(85)	-5.1(83)	-8.4(91)	-11.4(103)
11	100	-22.8(47)	-32.6(41)	-42.0(42)	-47.7(47)
	300	-14.7(65)	-20.0(77)	-29.6(68)	-34.3(75)
	500	-11.3(76)	-14.6(96)	-24.3(90)	-29.3(97)

6.2 Winsorization in a stratified simple random sampling without-replacement design

We also tested the calibration method described in Section 5. We generated a population of size $N=\text{5,000},$ which we divided into five strata, $U_1\mathrm{,}\dots \mathrm{,}{U}_5,$ of size $N_1\mathrm{,}\dots \mathrm{,}{N}_5\mathrm{,}$ respectively; see Table 6.4 for the values of $N_h.$ In each stratum, we generated a variable of interest $y$ according to a lognormal distribution with parameters $\text{log} \left(2,000\right)$ and $\mathrm{1.5.}$

From the population we selected $M=\text{5,000}$ samples according to a stratified simple random sampling without-replacement design. In stratum $U_h\mathrm{,}$ we selected a sample $S_h$ of size $n_h$ according to a simple random sampling without-replacement design; see Table 6.4 for the sizes $n_h$ and the corresponding sampling fractions, $f_h=n_h/N_h.$

The objective here is to estimate the total in the population, $t={\displaystyle {\sum }_{i\in U}}{y}_i\mathrm{,}$ and the stratum totals $t_h={\displaystyle {\sum }_{i\in U_h}{y}_i}\mathrm{,}h=1,\dots \mathrm{,}H\mathrm{.}$ In other words, in our example, the strata correspond to domains of interest. Since the strata form a partition of the population, we have the consistency relation, $t={\displaystyle {\sum }_{h=1}^H{t}_h}\mathrm{.}$ Similarly, the expansion estimators satisfy the consistency relation $\widehat{t}={\displaystyle {\sum }_{h=1}^H{\widehat{t}}_h}\mathrm{,}$ where $\widehat{t}={\displaystyle {\sum }_{i\in S}{d}_i{y}_i}$ and ${\widehat{t}}_h={\displaystyle {\sum }_{i\in S_h}{d}_i{y}_i}$ with $d_i=N_h/n_h$ if $i\in U_h\mathrm{.}$

For each sample, we first computed the robust estimator (4.8) in each stratum and aggregated the robust estimates to produce an aggregate robust estimate, ${\widehat{t}}_{R\left(\text{agg}\right)}={\displaystyle {\sum }_{h=1}^H{\widehat{t}}_{R\mathrm{,}h}}.$ Independently, we computed the robust estimator (4.8), denoted ${\widehat{t}}_{R\mathrm{,0}}\mathrm{,}$ at the population level. To ensure that the consistency relation (5.1) was satisfied, we performed the calibration procedure described in Section 5 to obtain the final robust estimates ${\widehat{t}}_{R\mathrm{,}h}^{*}\mathrm{,}h=0,\dots \mathrm{,5.}$ We used four systems of coefficients $q_h:$ (1) $q_0=0$ and $q_1=\dots =q_5=1;$ (2) $q_0=0$ and $q_h=n_h^{-1}\left(1-f_h\right)\mathrm{,}$ $h=1,\dots \mathrm{,5};$ (3) $q_0=0$ and $q_h=\text{CV}\left({\widehat{t}}_h\right)=$ $\sqrt{N_h^2\left(1-f_h\right)n_h^{-1}{S}_h^2}/t_h\mathrm{,}$ where $S_h^2={\left(N_h-1\right)}^{-1}{\displaystyle {\sum }_{i\in U_h}}{\left(y_i-{\overline{y}}_{Uh}\right)}^2\mathrm{,}h=1,\dots \mathrm{,5};$ (4) $q_0=0$ and $q_h=\widehat{\text{CV}}\left({\widehat{t}}_h\right)=\sqrt{N_h^2\left(1-f_h\right)n_h^{-1}{s}_h^2}/{\widehat{t}}_h\mathrm{,}$ where $s_h^2={\left(n_h-1\right)}^{-1}{\displaystyle {\sum }_{i\in S_h}}{\left(y_i-{\overline{y}}_{Sh}\right)}^2\mathrm{,}h=1,\dots \mathrm{,5.}$ We make the following remarks on the choice of the coefficients $q_h:$ (i) For all four systems, we assigned a weight $q_0=0$ to estimate ${\widehat{t}}_{R\mathrm{,0}},$ which is equivalent to making no change in the robust estimate at the population level. In other words, we have ${\widehat{t}}_{R\mathrm{,0}}^{*}={\widehat{t}}_{R\mathrm{,0}}.$ (ii) The first weighting system assigns an equal weight to all strata regardless of the sample size or sampling fraction. (iii) In the case of the second system, the coefficient $q_h$ is a function of the sample size $n_h$ and the sampling fraction $f_h,$ but it is independent of the intra-stratum variability $S_h^2.$ (iv) In the third and fourth systems, the choice of $q_h$ depends on the actual CV and the estimated CV respectively, for the reasons mentioned in Section 5.

Table 6.4
Characteristics of the strata
Table summary
This table displays the results of Characteristics of the strata. The information is grouped by Stratum (appearing as row headers), 1, 2, 3, 4 and 5 (appearing as column headers).

Stratum	1	2	3	4	5
$N_h$	2,000	1,500	1,000	400	100
$n_h$	20	75	100	80	80
$f_h$	0.01	0.05	0.1	0.2	0.8

For each robust estimator, we computed the Monte Carlo relative bias (as a percentage) and the relative efficiency (with respect to the expansion estimator); see Section 6.1. The results are presented in Table 6.5.

The results show that the initial robust estimators ${\widehat{t}}_{R\mathrm{,}h}$ are biased, as expected. The bias is larger in strata with a small sampling fraction. For example, in Stratum 1, for which $f_1=1\%,$ the relative bias of ${\widehat{t}}_{\mathrm{1,}h}$ is $-11.9\%,$ compared with only $-1.5\%$ in Stratum 5, for which $f_5=80\%.$ We also note that the initial robust estimators are all more efficient than the corresponding expansion estimator, with relative efficiency values ranging from 57% to 97%. The aggregate estimator ${\widehat{t}}_{R\left(\text{agg}\right)}$ obtained by summing the initial estimators ${\widehat{t}}_{R\mathrm{,}h}\mathrm{,}h=1,\dots \mathrm{,5}$ shows a modest bias with a value equal to $-5.7\%$ but is more efficient than the population-level expansion estimator $\widehat{t},$ with a relative efficiency of 87%.

The population-level winsorized estimator, ${\widehat{t}}_{R\mathrm{,0}},$ shows a small bias with a value equal to $-2.8\%$ and is significantly more efficient than the expansion estimator, with a relative efficiency of 81%. The final estimators ${\widehat{t}}_{R\mathrm{,}h}^{*}$ obtained using the system of coefficients $q_h=1$ for $h=1,\dots ,5$ all have lower bias than the initial estimator ${\widehat{t}}_{R\mathrm{,}h}\mathrm{,}$ except for Stratum 5. This is due to the fact that we force the sum of the final estimates ${\widehat{t}}_{R\mathrm{,}h}^{*}$ to calibrate on a low-bias estimator. On the other hand, the decrease in the bias is accompanied by a slight decrease in efficiency. For example, in Stratum 4, the relative efficiency is 63% for the robust estimator ${\widehat{t}}_{R\mathrm{,4}}$ and 66% for the final estimator ${\widehat{t}}_{R\mathrm{,4}}^{*}\mathrm{.}$ In the case of Stratum 5, the first system of coefficients is clearly unsuitable, since it leads to a change in the estimate for this stratum, like all the other strata, when this stratum has a very high sampling fraction of 80%. In fact, for this system of coefficients, the estimator ${\widehat{t}}_{R\mathrm{,5}}^{*}$ is less efficient than the expansion estimator, with a relative efficiency of 104. The second choice of coefficients $q_h,$ which takes the sampling fraction $f_h$ and the sample size $n_h$ into account, leads to some interesting results. The final robust estimator in Stratum 1, ${\widehat{t}}_{R\mathrm{,1}}^{*},$ has an appreciably lower bias than the initial estimator ${\widehat{t}}_{R\mathrm{,1}}$ and the final estimator based on the first system of coefficients, at the cost of a slight loss of efficiency. For Stratum 5, the estimator ${\widehat{t}}_{R\mathrm{,5}}^{*}$ has a low bias (a relative bias of $-0.8\%)$ and the same 97% efficiency as the initial estimator ${\widehat{t}}_{R\mathrm{,5}}.$ The third and fourth $q_h$ weighting systems lead to similar relative bias and relative efficiency results. For Stratum 1, they lead to lower relative biases than the first weighting system, at the cost of a slight loss of efficiency. For Strata 2, 3 and 4, all four systems of coefficients exhibit similar relative bias and relative efficiency. For Stratum 5, the final estimators are virtually unbiased and no less efficient that the expansion estimator.

Table 6.5
Monte Carol relative bias (in %) and relative efficiency (in parentheses) of the robust estimators at the global level and the stratum level
Table summary
This table displays the results of Monte Carol relative bias (in %) and relative efficiency (in parentheses) of the robust estimators at the global level and the stratum level. The information is grouped by Global estimator (appearing as row headers), ${\widehat{t}}_{R\left(\text{agg}\right)}$ , ${\widehat{t}}_{R\mathrm{,0}}={\widehat{t}}_{R\mathrm{,0}}^{*}$ , ${\widehat{t}}_{R\mathrm{,0}}={\widehat{t}}_{R\mathrm{,0}}^{*}$ , ${\widehat{t}}_{R\mathrm{,0}}={\widehat{t}}_{R\mathrm{,0}}^{*}$ and ${\widehat{t}}_{R\mathrm{,0}}={\widehat{t}}_{R\mathrm{,0}}^{*}$ (appearing as column headers).

Global estimator		${\widehat{t}}_{R\left(\text{agg}\right)}$	${\widehat{t}}_{R\mathrm{,0}}={\widehat{t}}_{R\mathrm{,0}}^{*}$	${\widehat{t}}_{R\mathrm{,0}}={\widehat{t}}_{R\mathrm{,0}}^{*}$	${\widehat{t}}_{R\mathrm{,0}}={\widehat{t}}_{R\mathrm{,0}}^{*}$	${\widehat{t}}_{R\mathrm{,0}}={\widehat{t}}_{R\mathrm{,0}}^{*}$
		-5.7(87)	-2.8(81)	-2.8(81)	-2.8(81)	-2.8(81)
		${\widehat{t}}_{R\mathrm{,}h}$	${\widehat{t}}_{R\mathrm{,}h}^{*}\left(q_h=1\right)$	${\widehat{t}}_{R\mathrm{,}h}^{*}\left(q_h=n_h^{-1}\left(1-f_h\right)\right)$	${\widehat{t}}_{R\mathrm{,}h}^{*}\left(q_h=\text{CV}\left({\widehat{t}}_h\right)\right)$	${\widehat{t}}_{R\mathrm{,}h}^{*}\left(q_h=\widehat{\text{CV}}\left({\widehat{t}}_h\right)\right)$
Stratum	1	-11.9(57)	-9.1(60)	-0.9(67)	-5.7(62)	-6.7(64)
	2	-6.3(74)	-3.4(76)	-3.3(76)	-3.3(76)	-3.1(78)
	3	-6.0(69)	-3.1(70)	-3.8(69)	-3.2(70)	-3.2(70)
	4	-6.6(63)	-3.7(66)	-4.2(65)	-3.3(66)	-3.4(70)
	5	-1.5(97)	1.5(104)	-0.8(97)	-0.2(98)	0.1(99)

Previous | Next