6. Simulation studies

Alina Matei and M. Giovanna Ranalli

Previous | Next

We evaluate the performance of the estimator presented in Section 5 by means of a Monte Carlo simulation under two different settings. The first one uses a real data set as the population and considers variables of interest that are all binary, while the second one uses simulated population data with variables of interest that are continuous. Results from the first setting are presented in Section 6.1, while those from the second setting are presented in Section 6.2.

In both settings, simple random sampling without replacement is employed and the following estimators are considered:

• $\text{HT=}{\displaystyle {\sum }_{k\in s}{y}_{kj}/{\pi }_k}:$ the Horvitz-Thompson estimator in the case of full response is computed as a benchmark in the absence of nonresponse.

• ${\widehat{Y}}_{j,\text{naive}}:$ the naive estimator given in (5.1); no explicit action is taken to adjust for unit and item nonresponse. Note that for simple random sampling without replacement, it reduces to ${\widehat{Y}}_{j,\text{naive}}=N{\displaystyle {\sum }_{k\in r_j}{y}_{kj}/n_{r_j}},$ where $n_{r_j}$ is the size of the set $r_j\mathrm{,}$ and it is the same as the Horvitz-Thompson estimator adjusted for unit nonresponse that assumes uniform response probabilities estimated by $n_{r_j}/n\mathrm{.}$

• ${\widehat{Y}}_{j,pq}:$ the three-phase estimator proposed in Section 5, Equation (5.2).

• ${\widehat{Y}}_{j,pq,\text{true}}:$ the three-phase estimator that uses the true values for the response probabilities $p_k$ and $q_{kj}$ is also computed for comparison with ${\widehat{Y}}_{j,pq}$ to understand the effect of estimating the response probabilities.

The simulations are carried out in R version 2.15, using the R package ‘ltm’ (Rizopoulos 2006) to fit the latent trait models. The following performance measures are computed for each estimator, generically denoted below by $\widehat{Y}$ where suffix $j$ is dropped for ease of notation $(Y$ denotes the population total):

• the Monte Carlo Bias

  $$\text{B=}{E}_{\text{sim}}\left(\widehat{Y}\right)-Y\mathrm{,}$$

where $E_{\text{sim}}\left(\widehat{Y}\right)={\displaystyle {\sum }_{i=1}^M{\widehat{Y}}_i}/M,{\widehat{Y}}_i$ is the value of the estimator $\widehat{Y}$ at the $i^{\text{th}}$ simulation run and $M$ is total number of simulation runs;

• the Relative Bias

  $$\text{RB=}\frac{\text{B}}{Y};$$

• the Monte Carlo Standard Deviation

  $$\sqrt{\text{VAR}}=\sqrt{\frac{1}{M-1}{\displaystyle \sum _{i=1}^M}{\left({\widehat{Y}}_i-E_{\text{sim}}\left(\widehat{Y}\right)\right)}^2};$$

• the Monte Carlo Mean Squared Error

  $${\text{MSE=B}}^2+\text{VAR}\mathrm{.}$$

6.1 Simulation setting 1

We consider the Abortion data set formed by four binary variables extracted from the 1986 British Social Attitudes Survey and concerning the attitude towards abortion. The data is available in the R package ‘ltm’ (Rizopoulos 2006). $N=\text{379}$ individuals answered the following questions after being asked if the law should allow abortion under the circumstances presented under each item:

1. The woman decides on her own that she does not wish to keep the baby.
2. The couple agrees that they do not wish to have a child.
3. The woman is not married and does not wish to marry the man.
4. The couple cannot afford any more children.

The variable of interest $y_j$ is selected to be the second one $\left(j=2\right)$ with a total $Y_j\text{=225}$ in the population.

The data is analyzed by Bartholomew et al. (2002) as an example in which a latent variable can be found that measures the attitude towards abortion. At the population level, we compute the latent variable (denoted here by ${\theta }_k^a)$ using Model (4.2) on the ${\left\{y_{k\ell }\right\}}_{k=1,\dots \mathrm{,}N;\ell =1,\dots \mathrm{,4}}$ data. The correlation between the values $y_{k\ell }$ and ${\theta }_k^a$ is approximatively equal to 0.85, for $\ell =1,\dots \mathrm{,4.}$ Afterwards, we have set ${\theta }_k={\widehat{\theta }}_k^a\mathrm{,}$ for all $k=1,\dots \mathrm{,}N\mathrm{.}$

At the population level, the unit response probabilities are generated using the following response model

$$p_k=1/\left(1+\exp \left(-\left(\text{0}\text{.7}+y_{k2}+{\theta }_k+\text{0}\text{.2}{\epsilon }_k\right)\right)\right)\mathrm{,}(6.1)$$

with ${\epsilon }_k\sim U\left(\mathrm{0,1}\right),$ to simulate nonignorable nonresponse. The population mean of $p_k$ is approximately 0.74.

To generate item response probabilities at the population level, the following model is used

$$q_{k\ell }=1/\left(1+\exp \left(-\left(b_{\ell }{\theta }_k+a_{\ell }+y_{k\ell }\right)\right)\right)\mathrm{,}\text{for }\ell =1,\dots \mathrm{,4,}(6.2)$$

where $b_{\ell }=3,$ for $\ell =1,\dots \mathrm{,4},$ while $a_{\ell }$ takes different values according to $\ell ;$ in particular, $a_1=1,a_2=0,a_3=-0.5$ and $a_4=1.$ The nominal item nonresponse rate for the four items in the population is 35%, 42%, 47%, 31%, respectively.

We draw $M=\text{10,000}$ simple random samples without replacement from the population using two sample sizes: $n=\text{50}$ and $n=\text{100}\mathrm{.}$ In each sample $s\mathrm{,}$ the units are classified as respondents according to Poisson sampling, using the probabilities $p_k$ computed as in Equation (6.1) and resulting in the set $r.$ Then, given $r,$ the matrix ${\left\{x_{k\ell }\right\}}_{k\in r;\ell =1,\mathrm{...},4}$ is constructed where the values $x_{k\ell }$ are drawn using Poisson sampling with probabilities $q_{k\ell }$ defined in (6.2). In each simulation run, Model (4.2) and the respondents set $r$ are used to compute the variable ${\widehat{\theta }}_k$ for all $k\in s$ as described in Section 4.4. Model (4.4) is fitted to obtain ${\widehat{p}}_k\mathrm{.}$ The average item nonresponse rate over simulations for the four items is found to be 26%, 33%, 38% and 23%. The jackknife variance estimator was computed as described in Section 5 using the gencalib() function in R package ‘sampling’ (Tillé and Matei 2012) and the logistic distance (Deville, Särndal and Sautory 1993).

Table 6.1 reports the results for $n=\text{50}$ and $n=\text{100}.$ As expected, $\text{HT}$ and ${\widehat{Y}}_{j,pq,\text{true}}$ have almost zero bias, with the second one showing a relatively larger MSE that is due uniquely to the smaller sample size. The naive estimator shows a very large negative bias. This is due to the fact that units with a zero value of $y_j$ are less likely to respond and the total is clearly underestimated. The estimator ${\widehat{Y}}_{j,pq}$ shows a much smaller bias than the naive estimator. Note that the performance of the proposed estimator is mostly driven by absolute bias, so that the performance is not particularly different when increasing the sample size, apart from a decrease in variance. If we compare ${\widehat{Y}}_{j,pq,\text{true}}$ and ${\widehat{Y}}_{j,pq},$ we note that ${\widehat{Y}}_{j,pq}$ still suffers from some bias that comes from response model misspecification (we are not accounting for the variables of interest values).

For the proposed estimator, the jackknife variance estimator was also tested by looking at the empirical coverage of a 95% confidence interval computed for each replicate as ${\widehat{Y}}_{j,pq}\pm \text{1}\text{.96}\sqrt{{\widehat{V}}_r}.$ For $n=\text{50}\mathrm{,}$ the mean value of $\sqrt{{\widehat{V}}_r}$ over simulations was 54.8, while for $n=\text{100}\mathrm{,}$ 53.3, with a 95% coverage rate of 94.6% and 96.3%, respectively. The replicate estimator overestimates the Monte Carlo standard deviation reported for ${\widehat{Y}}_{j,pq}$ in Table 6.1 in both cases, but shows good coverage rates.

Table 6.1
Simulation results for setting 1 - Abortion data set
Table summary
This table displays the results of Simulation results for setting 1 - Abortion data set. The information is grouped by Estimator (appearing as row headers), B, $\sqrt{\text{VAR}}$ , MSE and % RB (appearing as column headers).

Estimator	B	$\sqrt{\text{VAR}}$	MSE	% RB
$n=\text{50}$
$\text{HT}$	0.05	24.5	600.5	$<$ 0.1
${\widehat{Y}}_{j,\text{naive}}$	-126.5	19.4	16,378.6	-56.2
${\widehat{Y}}_{j,pq}$	20.6	32.4	1,474.1	9.1
${\widehat{Y}}_{j,pq,\text{true}}$	0.02	35.0	1,225.0	$<$ 0.1
$n=\text{100}$
$\text{HT}$	-0.06	16.0	255.5	$<$ 0.1
${\widehat{Y}}_{j,\text{naive}}$	-126.9	13.5	16,284.1	-56.4
${\widehat{Y}}_{j,pq}$	17.9	21.9	802.2	8.0
${\widehat{Y}}_{j,pq,\text{true}}$	-0.1	23.7	559.9	$<$ 0.1

To study the performance of the latent model on the population level and the correlation between the variable of interest and the estimated latent variable, we apply the procedure described earlier using $q_{k\ell }$ defined in (6.2) to construct the matrix ${\left\{x_{k\ell }\right\}}_{k=1,\dots \mathrm{,}N;\ell =1,\mathrm{...},4}$ for all population units. We fit Model (4.2) on the population level and compute the variable ${\theta }_k$ for all $k=1,\dots \mathrm{,}N.$ The Cronbach’s alpha measure takes value 0.83 showing a good internal consistency of the items. The correlation coefficient between the variable of interest and the estimated latent variable takes value 0.76, indicating that the latent auxiliary information has a strong power of predicting $y_{k2}\mathrm{,}$ as advocated in the model of Cassel et al. (1983). Inspection of the two-way margins for the matrix $\left\{x_{k\ell }\right\}$ gives the residuals ${\left(O-E\right)}^2/E$ between 0.03 and 0.23. Similarly, the three-way margins for the matrix $\left\{x_{k\ell }\right\}$ give residuals between 0 and 1.19. This indicates that we have no reason to reject here the one-factor latent Model (4.2) (see Bartholomew et al. 2002, page 186).

6.2 Simulation setting 2

We generate $\left\{y_{k1}\mathrm{,}\dots \mathrm{,}{y}_{k6}\mathrm{,}{\theta }_k\right\}$ for $k=1,\dots \mathrm{,}N=\text{2,000}$ using a multivariate normal distribution with mean 1. The degree of correlation between $y_{\ell }$ and $y_{{\ell }^{\prime }}$ is 0.8, with $\ell \mathrm{,}{\ell }^{\prime }=1,\dots \mathrm{,6,}\ell \ne {\ell }^{\prime }.$ We set the variable of interest to be $y_6$ and consider different degrees of correlation between its values and those taken by ${\theta }_k,$ namely 0.3, 0.5, 0.8. The values of ${\theta }_k$ are afterwards standardized to have mean 0 and variance 1.

The response probabilities are obtained by first computing

$$p_k^{\circ }=1/\left[1+\exp \left(-\left(0.5+y_{k1}+{\theta }_k\right)\right)\right]\mathrm{,}\text{for }k=1,\dots \mathrm{,}N\mathrm{,}(6.3)$$

and then rescaling them to take values between 0.1 and 0.9, with a population mean approximatively equal to 0.7.

The item response probabilities are generated by first computing

$$q_{k\ell }^{\circ }=1/\left(1+\exp \left(-\left(b_{\ell }{\theta }_k+a_{\ell }+y_{k\ell }\right)\right)\right)\mathrm{,}\text{for }k=1,\dots \mathrm{,}N\text{ and }\ell =1,\dots \mathrm{,6,}(6.4)$$

where ${\left\{a_{\ell }\right\}}_{\ell =1,\dots \mathrm{,6}}=\left\{\mathrm{1,0,}-\text{0}\text{.5}\mathrm{,1,0,}-\text{0}\text{.5}\right\}$ and ${\left\{b_{\ell }\right\}}_{\ell =1,\dots \mathrm{,6}}=\left\{\mathrm{1,1,1,}\text{1}\text{.5}\mathrm{,}\text{1}\text{.5}\mathrm{,}\text{1}\text{.5}\right\}\mathrm{,}$ and then rescaling the values to be between 0.1 and 0.95.

We draw $M=\text{10,000}$ samples by simple random sampling without replacement of size $n=\text{200}.$ For each sample $s,$ a response set $r$ is created by carrying out Poisson sampling with parameter $p_k$ defined in (6.3). Each element of the matrix ${\left\{x_{k\ell }\right\}}_{k\in r\mathrm{,}\ell =1,\dots \mathrm{,6}}$ is generated using Poisson sampling with parameter $q_{k\ell }$ defined in (6.4). Item nonresponse rates over simulations take approximately value 18%, 28%, 35%, 19%, 29%, 34%, for $\ell =1,\dots \mathrm{,6},$ respectively. For each simulation run, Model (4.2) is used to compute the variable ${\widehat{\theta }}_k$ for all $k\in s.$ Model (4.4) is then fitted to obtain ${\widehat{p}}_k\mathrm{.}$

Table 6.2
Simulation results for setting 2 - Simulated continuous data
Table summary
This table displays the results of Simulation results for setting 2 - Simulated continuous data. The information is grouped by Estimator (appearing as row headers), B, $\sqrt{\text{VAR}}$, MSE and % RB (appearing as column headers).

Estimator	B	$\sqrt{\text{VAR}}$	MSE	% RB
correlation coefficient 0.3
$\text{HT}$	-0.7	131.6	17,331.2	$\approx$ -0.0
${\widehat{Y}}_{j\mathrm{,}\text{naive}}$	825.6	177.1	713,039.3	41.0
${\widehat{Y}}_{j\mathrm{,}pq}$	-227.4	188.0	87,033.0	-11.3
${\widehat{Y}}_{j\mathrm{,}pq\mathrm{,}\text{true}}$	48.4	231.8	56,073.2	2.4
correlation coefficient 0.5
$\text{HT}$	0.1	135.0	18,220.5	$\approx$ 0.0
${\widehat{Y}}_{j\mathrm{,}\text{naive}}$	972.6	176.2	977,009.5	50.7
${\widehat{Y}}_{j\mathrm{,}pq}$	-180.0	175.5	63,552.0	-9.4
${\widehat{Y}}_{j\mathrm{,}pq\mathrm{,}\text{true}}$	74.8	212.7	50,844.0	3.9
correlation coefficient 0.8
$\text{HT}$	-0.1	134.1	17,992.0	$\approx$ -0.0
${\widehat{Y}}_{j\mathrm{,}\text{naive}}$	1,154.6	168.1	1,361,388.1	57.7
${\widehat{Y}}_{j\mathrm{,}pq}$	-184.8	164.4	61,173.0	-9.2
${\widehat{Y}}_{j\mathrm{,}pq\mathrm{,}\text{true}}$	100.6	196.2	48,597.9	5.0

Table 6.2 reports on the performance of the estimators for the three values taken by the nominal correlation coefficient between $y_{k1}$ and ${\theta }_k\text{:0}\text{.3},\text{0}\text{.5}$ and $\text{0}\text{.8}.$ The proposed estimator is always able to reduce bias over the naive estimator, even when the correlation between the variable of interest and the latent variable gets smaller. The relative bias takes acceptable values in most cases. Bias deserves a closer look. The naive estimator in all cases largely overestimates the total. This is expected, because the values $p_k,q_{k6},{\theta }_k$ and $y_{k6}$ all go in the same direction. Therefore, in our respondents sample, we are more likely to find relative larger values for $y_6$ by this providing overestimation for the naive estimator. On the other hand, ${\widehat{Y}}_{j,pq}$ underestimates the total because it is based only on the observed units of $r_j$ that do have relatively large values for $y_6,$ but also relatively large values for $p_k$ and $q_{k6}$ and, therefore, end up having a small weight.

The matrix of population values ${\left\{x_{k\ell }\right\}}_{k=1;\dots ;\text{2,000;}\ell =1,\dots \mathrm{,6}}$ is constructed in the same way as in Section 6.1 to validate the assumptions behind the 2PL model. The Cronbach’s alpha takes approximately value 0.5 for the correlation coefficient equal to 0.3, 0.6 for 0.5, and 0.7 for 0.8; the pairwise association between the six items reveals $p\text{-}$values smaller than 0.01. Inspection of the two-way and three-way margins of the matrix $\left\{x_{k\ell }\right\}$ gives residuals ${\left(O-E\right)}^2/E$ that all take values smaller than 4. Therefore, the one factor latent model can be accepted and items all seem to be measuring the same latent trait.

Previous | Next