Skip to contents

In the following we present the methodology in surveysd by applying the workflow described in vignette("surveysd") to multiple consecutive years of EU-SILC data for one country. The methodology contains the following steps, in this order

  • Draw BB bootstrap replicates from EU-SILC data for each year yty_t, t=1,,nyt=1,\ldots,n_y separately. Since EU-SILC has a rotating panel design the bootstrap replicate of a household is carried forward through the years. That is, the bootstrap replicate of a household in the follow-up years is set equal to the bootstrap replicate of the same household when it first enters EU-SILC.
  • Multiply each set of bootstrap replicates by the sampling weights to obtain uncalibrated bootstrap weights and calibrate each of the uncalibrated bootstrap weights using iterative proportional fitting.
  • Estimate the point estimate of interest θ\theta, for each year and each calibrated bootstrap weight to obtain θ̃(i,yt)\tilde{\theta}^{(i,y_t)}, t=1,,nyt=1,\ldots,n_y, i=1,,Bi=1,\ldots,B. For fixed yty_t apply a filter with equal weights for each ii on θ̃(i,y*)\tilde{\theta}^{(i,y^*)}, y*{yt1,yt,yt+1}y^*\in \{y_{t-1},y_{t},y_{t+1}\} , to obtain θ̃(i,yt)\tilde{\theta}^{(i,y_t)}. Estimate the variance of θ\theta using the distribution of θ̃(i,yt)\tilde{\theta}^{(i,y_t)}.

Bootstrapping

Bootstrapping has long been around and used widely to estimate confidence intervals and standard errors of point estimates.[Efron (1979)} Given a random sample (X1,,Xn)(X_1,\ldots,X_n) drawn from an unknown distribution FF the distribution of a point estimate θ(X1,,Xn;F)\theta(X_1,\ldots,X_n;F) can in many cases not be determined analytically. However when using bootstrapping one can simulate the distribution of θ\theta.

Let s(.)s_{(.)} be a bootstrap sample, e.g. drawing nn observations with replacement from the sample (X1,,Xn)(X_1,\ldots,X_n), then one can estimate the standard deviation of θ\theta using BB bootstrap samples through sd(θ)=1B1i=1B(θ(si)θ¯)2,sd(\theta) = \sqrt{\frac{1}{B-1}\sum\limits_{i=1}^B (\theta(s_i)-\overline{\theta})^2} \quad,

with θ¯:=1Bi=1Bθ(si)\overline{\theta}:=\frac{1}{B}\sum\limits_{i=1}^B\theta(s_i) as the sample mean over all bootstrap samples.

In context of sample surveys with sampling weights one can use bootstrapping to calculate so called bootstrap weights. These are computed via the bootstrap samples sis_{i}, i=1,,Bi=1,\ldots,B, where for each sis_{i} every unit of the original sample can appear 00- to mm-times. With fjif_j^{i} as the frequency of occurrence of observation jj in bootstrap sample sis_i the uncalibrated bootstrap weights b̃ji\tilde{b}_{j}^{i} are defined as:

b̃ji=fjiwj, \tilde{b}_{j}^{i} = f_j^{i} w_j \quad,

with wjw_j as the calibrated sampling weight of the original sample. Using iterative proportional fitting procedures one can recalibrate the bootstrap weights b̃j.\tilde{b}_{j}^{.}, j=1,,Bj=1,\ldots,B to get the adapted or calibrated bootstrap weights bjib_j^i, j=1,,Bj=1,\ldots,B.

Rescaled Bootstrap

Since EU-SILC is a stratified sample without replacement drawn from a finite population the naive bootstrap procedure, as described above, does not take into account the heterogeneous inclusion probabilities of each sample unit. Thus it will not yield satisfactory results. Therefore we will use the so called rescaled bootstrap procedure introduced and investigated by (Rao and Wu 1988). The bootstrap samples are selected without replacement and do incorporate the stratification as well as clustering on multiple stages (see (Chipperfield and Preston 2007),(Preston 2009)).

For simplistic reasons we will only describe the rescaled bootstrap procedure for a two stage stratified sampling design. For more details on a general formulation please see (Preston 2009).

Sampling design

Consider the finite population UU which is divided into HH non-overlapping strata h=1,,HUh=U\bigcup\limits_{h=1,\ldots,H} U_h = U, of which each strata hh contains of NhN_h clusters. For each strata hh, ChcC_{hc}, c=1,,nhc=1,\ldots,n_h clusters are drawn, containing NhcN_{hc} households. Furthermore in each cluster ChcC_{hc} of each strata hh simple random sampling is performed to select a set of households YhcjY_{hcj}, j=1,,nhcj=1,\ldots,n_{hc}.

Bootstrap procedure

In contrast to the naive bootstrap procedure where for a stage, containing nn sampling units, the bootstrap replicate is obtained by drawing nn sampling units with replacement, for the rescaled bootstrap procedure n*=n2n^*=\left\lfloor\frac{n}{2}\right\rfloor sampling units are drawn without replacement. Given a value xx, x\lfloor x\rfloor denotes the largest integer smaller than xx, whereas x\lceil x\rceil denotes the smallest integer lager then xx. (Chipperfield and Preston 2007) have shown that the choice of either n2\left\lfloor\frac{n}{2}\right\rfloor or n2\left\lceil\frac{n}{2}\right\rceil is optimal for bootstrap samples without replacement, although n2\left\lfloor\frac{n}{2}\right\rfloor has the desirable property that the resulting uncalibrated bootstrap weights will never be negative.

At the first stage the ii-th bootstrap replicate, fhci,1f^{i,1}_{hc}, for each cluster ChcC_{hc},c=1,,nhc=1,\ldots,n_h, belonging to strata hh, is defined by

fhci,1=1λh+λhnhnh*δhcc{1,,nh} f^{i,1}_{hc} = 1-\lambda_h+\lambda_h\frac{n_h}{n_h^*}\delta_{hc} \quad\quad \forall c \in \{1,\ldots,n_h\} with nh*=nh2 n_h^* = \left\lfloor\frac{n_h}{2}\right\rfloor λh=nh*(1nhNh)nhnh*, \lambda_h = \sqrt{\frac{n_h^*(1-\frac{n_h}{N_h})}{n_h-n_h^*}} \quad ,

where δhc=1\delta_{hc}=1 if cluster cc is selected in the sub-sample of size nh*n_h^* and 0 otherwise.

The ii-th bootstrap replicate at the second stage, fhcji,2f^{i,2}_{hcj}, for each household YhcjY_{hcj}, j=1,,nhcj=1,\ldots,n_{hc}, belonging to cluster cc in strata hh is defined by

fhcji,2=fhci,1λhcnhnh*δhc[nhcnhc*δhcj1]c{1,,nh} f^{i,2}_{hcj} = f^{i,1}_{hc} - \lambda_{hc}\sqrt{\frac{n_h}{n_h^*}}\delta_{hc}\left[\frac{n_{hc}}{n_{hc}^*}\delta_{hcj}-1\right] \quad\quad \forall c \in \{1,\ldots,n_h\} with nhc*=nhc2 n_{hc}^* = \left\lfloor\frac{n_{hc}}{2}\right\rfloor λhc=nhc*Nh(1nhcNhc)nhcnhc*, \lambda_{hc} = \sqrt{\frac{n_{hc}^*N_h(1-\frac{n_{hc}}{N_{hc}})}{n_{hc}-n_{hc}^*}} \quad ,

where δhcj=1\delta_{hcj}=1 if household jj is selected in the sub sample of size nhc*n_{hc}^* and 0 otherwise.

Single PSUs

When dealing with multistage sampling designs the issue of single PSUs, e.g. a single response unit is present at a stage or in a strata, can occur. When applying bootstrapping procedures these single PSUs can lead to a variety of issues. For the methodology proposed in this work we combined single PSUs at each stage with the next smallest strata or cluster, before applying the bootstrap procedure.

Taking bootstrap replicates forward

The bootstrap procedure above is applied on the EU-SILC data for each year yty_t, t=1,,nyt=1,\ldots,n_y separately. Since EU-SILC is a yearly survey with rotating penal design the ii-th bootstrap replicate at the second stage, fhcji,2f^{i,2}_{hcj}, for a household YhcjY_{hcj} is taken forward until the household YhcjY_{hcj} drops out of the sample. That is, for the household YhcjY_{hcj}, which enters EU-SILC at year y1y_1 and drops out at year yt̃y_{\tilde{t}}, the bootstrap replicates for the years y2,,yt̃y_2,\ldots,y_{\tilde{t}} are set to the bootstrap replicate of the year y1y_1.

Split households

Due to the rotating penal design so called split households can occur. For a household participating in the EU-SILC survey it is possible that one or more residents move to a new so called split household, which is followed up on in the next wave. To take this dynamic into account we extended the procedure of taking forward the bootstrap replicate of a household for consecutive waves of EU-SILC by taking forward the bootstrap replicate to the split household. That means, that also any new individuals in the split household will inherit this bootstrap replicate.

Taking bootstrap replicates forward as well as considering split households ensures that bootstrap replicates are more comparable in structure with the actual design of EU-SILC.

Uncalibrated bootstrap weights

Using the ii-th bootstrap replicates at the second stage one can calculate the ii-th uncalibrated bootstrap weights bhcjib_{hcj}^{i} for each household YhcjY_{hcj} in cluster cc contained in strata hh by

b̃hcji=fhcji,2whcj, \tilde{b}_{hcj}^{i} = f^{i,2}_{hcj} w_{hcj} \quad, where whcjw_{hcj} corresponds to the original household weight contained in the sample.

For ease of readability we will drop the subindices regarding strata hh and cluster cc for the following sections, meaning that the jj-th household in cluster cc contained in strata hh, YhcjY_{hcj}, will now be denoted as the jj-th household, YjY_{j}, where jj is the position of the household in the data. In accordance to this the ii-th uncalibrated bootstrap replicates for household jj are thus denoted as b̃ji\tilde{b}_j^{i} and the original household weight as wjw_j.

Iterative proportional fitting (IPF)

The uncalibrated bootstrap weights b̃ji\tilde{b}_j^{i} computed through the rescaled bootstrap procedure yields population statistics that differ from the known population margins of specified sociodemographic variables for which the base weights wjw_j have been calibrated. To adjust for this the bootstrap weights b̃ji\tilde{b}_{j}^{i} can be recalibrated using iterative proportional fitting as described in (Meraner, Gumprecht, and Kowarik 2016).

Let the original weight wjw_{j} be calibrated for n=nP+nHn=n_P+n_H sociodemographic variables which are divided into the sets 𝒫:={pc,c=1,nP}\mathcal{P}:=\{p_{c}, c=1 \ldots,n_P\} and :={hc,c=1,nH}\mathcal{H}:=\{h_{c}, c=1 \ldots,n_H\}. 𝒫\mathcal{P} and \mathcal{H} correspond to personal, for example gender or age, or household variables, like region or households size, respectively. Each variable in either 𝒫\mathcal{P} or \mathcal{H} can take on PcP_{c} or HcH_{c} values with and NvpcN^{p_c}_v, v=1,,Pcv=1,\ldots,P_c, or NvhcN^{h_c}_v, v=1,,Hcv=1,\ldots,H_c, as the corresponding population margins. Starting with k=0k=0 the iterative proportional fitting procedure is applied on each b̃ji\tilde{b}_j^{i}, i=1,,Bi=1,\ldots, B separately. The weights are first updated for personal and afterwards updated for household variables. If constraints regarding the populations margins are not met kk is raised by 1 and the procedure starts from the beginning. For the following denote as starting weight b̃j[0]:=b̃ji\tilde{b}_j^{[0]}:=\tilde{b}_j^{i} for fixed ii.

Adjustment and trimming for 𝒫\mathcal{P}

The uncalibrated bootstrap weight b̃j[(n+1)k+c1]\tilde{b}_j^{[(n+1)k+c-1]} for the jj-th observation is iteratively multiplied by a factor so that the projected distribution of the population matches the respective calibration specification NpcN_{p_c}, c=1,,nPc=1, \ldots,n_P. For each c{1,,nP}c \in \left\{1, \ldots,n_P\right\} the calibrated weights against NvpcN^{p_c}_v are computed as b̃j[(n+1)k+c]=b̃j[(n+1)k+c1]Nvpclb̃l[(n+1)k+c1], \tilde{b}_j^{[(n+1)k+c]} = {\tilde{b}_j}^{[(n+1)k+c-1]}\frac{N^{p_c}_v}{{\sum\limits_l} {\tilde{b}}_l^{[(n+1)k+c-1]}}, where the summation in the denominator expands over all observations which have the same value as observation jj for the sociodemographic variable pcp_c. If any weights b̃j[nk+c]\tilde{b}_j^{[nk+c]} fall outside the range [wj4;4wj]\left[\frac{w_j}{4};4w_j\right] they will be recoded to the nearest of the two boundaries. The choice of the boundaries results from expert-based opinions and restricts the variance of which has a positive effect on the sampling error. This procedure represents a common form of weight trimming where very large or small weights are trimmed in order to reduce variance in exchange for a possible increase in bias ((Potter 1990),(Potter 1993)).

Averaging weights within households

Since the sociodemographic variables p1,,pncp_1,\ldots,p_{n_c} include person-specific variables, the weights b̃j[nk+np]\tilde{b}_j^{[nk+n_p]} resulting from the iterative multiplication can be unequal for members of the same household. This can lead to inconsistencies between results projected with household and person weights. To avoid such inconsistencies each household member is assigned the mean of the household weights. That is for each person jj in household aa with hah_a household members, the weights are defined by b̃j[(n+1)k+np+1]=lab̃l[(n+1)k+np]ha \tilde{b}_j^{[(n+1)k+n_p+1]} = \frac{{\sum\limits_{l\in a}} {\tilde{b}_l^{[(n+1)k+n_p]}}}{h_a} This can result in losing the population structure performed in the previous subsection.

Adjustment and trimming for \mathcal{H}

After adjustment for individual variables the weights bj[nk+np+1]b_j^{[nk+n_p+1]} are updated for the set of household variables \mathcal{H} according to a household convergence constraint parameter ϵh\epsilon_h. The parameters ϵh\epsilon_h represent the allowed deviation from the population margins using the weights bj[nk+np+1]b_j^{[nk+n_p+1]} compared to NvhcN^{h_c}_v, c=1,,nHc=1,\ldots,n_H, v=1,,Hcv=1,\ldots,H_c. The updated weights are computed as bj[(n+1)k+np+c+1]={bj[(n+1)k+np+1]Nvhclbl[(n+1)k+np+1]if lbj[(n+1)k+np+1]((10.9ϵh)Nvhc,(1+0.9ϵh)Nvhc)bj[(n+1)k+np+1]otherwise b_j^{[(n+1)k+n_p+c+1]} = \begin{cases} b_j^{[(n+1)k+n_p+1]}\frac{N^{h_c}_v}{\sum\limits_{l} b_l^{[(n+1)k+n_p+1]}} \quad \text{if } \sum\limits_{l} b_j^{[(n+1)k+n_p+1]} \notin ((1-0.9\epsilon_h)N^{h_c}_v,(1+0.9\epsilon_h)N^{h_c}_v) \\ b_j^{[(n+1)k+n_p+1]} \quad \text{otherwise} \end{cases} with the summation in the denominator ranging over all households ll which take on the same values for hch_c as observation jj. As described in the previous subsection the new weight are recoded if they exceed the interval [wj4;4wj][\frac{w_j}{4};4w_j] and set to the upper or lower bound, depending of bj[(n+1)k+np+c+1]b_j^{[(n+1)k+n_p+c+1]} falls below or above the interval respectively.

Convergence

For each adjustment and trimming step the factor Nv(.)lbl[(n+1)k+j]\frac{N^{(.)}_v}{\sum\limits_{l} b_l^{[(n+1)k+j]}}, j{1,,n+1}{np+1}j\in \{1,\ldots,n+1\}\backslash \{n_p+1\}, is checked against convergence constraints for households, ϵh\epsilon_h, or personal variables ϵp\epsilon_p, where (.)(.) corresponds to either a household or personal variable. To be more precise for variables in 𝒫\mathcal{P} the constraints

Nvpclb̃l[(n+1)k+j]((1ϵp)Nvpc,(1+ϵp)Nvpc) \frac{N^{p_c}_v}{{\sum\limits_l} {\tilde{b}}_l^{[(n+1)k+j]}} \in ((1-\epsilon_p)N^{p_c}_v,(1+\epsilon_p)N^{p_c}_v) and for variables in \mathcal{H} the constraints

Nvhclb̃l[(n+1)k+j]((1ϵh)Nvhc,(1+ϵh)Nvhc) \frac{N^{h_c}_v}{{\sum\limits_l} {\tilde{b}}_l^{[(n+1)k+j]}} \in ((1-\epsilon_h)N^{h_c}_v,(1+\epsilon_h)N^{h_c}_v) are verified, where the sum in the denominator expands over all observations which have the same value for variables hch_c or pcp_c. If these constraints hold true the algorithm reaches convergence, otherwise kk is raised by 1 and the procedure repeats itself.

The above described calibration procedure is applied on each year yty_t of EU-SILC separately, t=1,nyt=1,\ldots n_y, thus resulting in so called calibrated bootstrap sample weights bj(i,yt)b_{j}^{(i,{y_t})}, i=1,,Bi=1,\ldots,B for each year yy and each household jj.

Variance estimation

Applying the previously described algorithms to EU-SILC data for multiple consecutive years yty_t, t=1,nyt=1,\ldots n_y, yields calibrated bootstrap sample weights bj(i,yt)b_{j}^{(i,{y_t})} for each year yty_t. Using the calibrated bootstrap sample weights it is straight forward to compute the standard error of a point estimate θ(𝐗yt,𝐰yt)\theta(\textbf{X}^{y_t},\textbf{w}^{y_t}) for year yty_t with 𝐗yt=(X1yt,,Xnyt)\textbf{X}^{y_t}=(X_1^{y_t},\ldots,X_n^{y_t}) as the vector of observations for the variable of interest in the survey and 𝐰yt=(w1yt,,wnyt\textbf{w}^{y_t}=(w_1^{y_t},\ldots,w_n^{y_t} as the corresponding weight vector, with

sd(θ)=1B1i=1B(θ(i,yt)θ(.,yt)¯)2 sd(\theta) = \sqrt{\frac{1}{B-1}\sum\limits_{i=1}^B (\theta^{(i,y_t)}-\overline{\theta^{(.,y_t)}})^2} with θ(.,yt)¯=1Bi=1Bθ(i,yt), \overline{\theta^{(.,y_t)}} = \frac{1}{B}\sum\limits_{i=1}^B\theta^{(i,y_t)} \quad, where θ(i,yt):=θ(𝐗yt,𝐛(i,yt))\theta^{(i,y_t)}:=\theta(\textbf{X}^{y_t},\textbf{b}^{(i,{y_t})}) is the estimate of θ\theta in the year yty_t using the ii-th vector of calibrated bootstrap weights.

As already mentioned the standard error estimation for indicators in EU-SILC yields high quality results for NUTS1 or country level. When estimation indicators on regional or other sub-aggregate levels one is confronted with point estimates yielding high variance.

To overcome this issue we propose to estimate θ\theta for 3, consecutive years using the calibrated bootstrap weights, thus calculating {θ(i,yt1),θ(i,yt),θ(i,yt+1)}\{\theta^{(i,y_{t-1})},\theta^{(i,y_t)},\theta^{(i,y_{t+1})}\}, i=1,,Bi=1,\ldots,B. For fixed ii one can apply a filter with equal filter weights on the time series {θ(i,yt1),θ(i,yt),θ(i,yt+1)}\{\theta^{(i,y_{t-1})},\theta^{(i,y_t)},\theta^{(i,y_{t+1})}\} to create θ̃(i,yt)\tilde{\theta}^{(i,y_t)}

θ̃(i,yt)=13[θ(i,yt1)+θ(i,yt)+θ(i,yt+1)]. \tilde{\theta}^{(i,y_t)} = \frac{1}{3}\left[\theta^{(i,y_{t-1})}+\theta^{(i,y_t)}+\theta^{(i,y_{t+1})}\right] \quad .

Doing this for all ii, i=1,,Bi=1,\ldots,B, yields θ̃(i,yt)\tilde{\theta}^{(i,y_t)}, i=1,,Bi=1,\ldots,B. The standard error of θ\theta can then be estimated with

sd(θ)=1B1i=1B(θ̃(i,yt)θ̃(.,yt)¯)2 sd(\theta) = \sqrt{\frac{1}{B-1}\sum\limits_{i=1}^B (\tilde{\theta}^{(i,y_t)}-\overline{\tilde{\theta}^{(.,y_t)}})^2} with θ̃(.,yt)¯=1Bi=1Bθ̃(i,yt). \overline{\tilde{\theta}^{(.,y_t)}}=\frac{1}{B}\sum\limits_{i=1}^B\tilde{\theta}^{(i,y_t)} \quad.

Applying the filter over the time series of estimated θ(i,yt)\theta^{(i,y_t)} leads to a reduction of variance for θ\theta since the filter reduces the noise in {θ(i,yt1),θ(i,yt),θ(i,yt+1)}\{\theta^{(i,y_{t-1})},\theta^{(i,y_t)},\theta^{(i,y_{t+1})}\} and thus leading to a more narrow distribution for θ̃(i,yt)\tilde{\theta}^{(i,y_t)}.

It should also be noted that estimating indicators from a survey with rotating panel design is in general not straight forward because of the high correlation between consecutive years. However with our approach to use bootstrap weights, which are independent from each other, we can bypass the cumbersome calculation of various correlations, and apply them directly to estimate the standard error. (Bauer et al. 2013) showed that using the proposed method on EU-SILC data for Austria the reduction in resulting standard errors corresponds in a theoretical increase in sample size by about 25%\%. Furthermore this study compared this method to the use of small area estimation techniques and on average the use of bootstrap sample weights yielded more stable results.

References

Bauer, Martin, Matthias Till, Richard Heuberger, Marcel Bilgili, Thomas Glaser, Elisabeth Kafka, Johannes Klotz, et al. 2013. “Studie Zu Armut Und Sozialer Eingliederung in Den Bundesl"andern.” Statistik Austria [in German].
Chipperfield, James, and John Preston. 2007. “Efficient Bootstrap for Business Surveys.” Survey Methodology 33 (December): 167–72. https://www150.statcan.gc.ca/n1/en/catalogue/12-001-X200700210494.
Efron, B. 1979. “Bootstrap Methods: Another Look at the Jackknife.” Ann. Statist. 7 (1): 1–26. https://doi.org/10.1214/aos/1176344552.
Meraner, Angelika, Daniela Gumprecht, and Alexander Kowarik. 2016. “Weighting Procedure of the Austrian Microcensus Using Administrative Data.” Austrian Journal of Statistics 45 (June): 3. https://doi.org/10.17713/ajs.v45i3.120.
Potter, Frank J. 1990. “A Study of Procedures to Identify and Trim Extreme Sampling Weights.” Proceedings of the American Statistical Association, Section on Survey Research Methods, 225–30. http://www.asasrms.org/Proceedings/papers/1990_034.pdf.
———. 1993. “The Effect of Weight Trimming on Nonlinear Survey Estimates.” Proceedings of the American Statistical Association, Section on Survey Research Methods 2: 758–63. http://www.asasrms.org/Proceedings/papers/1993_127.pdf.
Preston, J. 2009. “Rescaled Bootstrap for Stratified Multistage Sampling.” Survey Methodology 35 (December): 227–34. https://www150.statcan.gc.ca/n1/en/catalogue/12-001-X200900211044.
Rao, J. N. K., and C. F. J. Wu. 1988. “Resampling Inference with Complex Survey Data.” Journal of the American Statistical Association 83 (401): 231–41.