Abstract
Surveys with a rotating panel design are a prominent tool for producing more efficient estimates for indicators regarding trends or net changes over time. Variance estimation for net changes becomes however more complicated due to a possibly high correlation between the panel waves. Therefore, these estimates are quite burdensome to produce with traditional means. With the R-package surveysd, we present a tool which supports a straightforward way for producing estimates and corresponding standard errors for complex surveys with a rotating panel design. The package uses bootstrap techniques which incorporate the panel design and thus makes it easy to estimate standard errors. In addition the package supports a method for producing more efficient estimates by cumulating multiple consecutive sample waves. This method can lead to a significant decrease in variance assuming that structural patterns for the indicator in question remain fairly robust over time. The usability of the package and variance improvement, using this bootstrap methodology, is demonstrated on data from the user database (UDB) for the EU Statistics on Income and Living Conditions of selected countries with various sampling designs.
Introduction
Many sociodemographic population parameters are periodically estimated, for instance the unemployment rate, wealth distribution or poverty rate. One important aspect for the regular estimation of those parameters is to observe net changes over time and to consequently monitor the effect of related policies. Surveys with a rotating panel design are a prominent tool for measuring net changes of an indicator over time. They are a compromise between a single panel survey, which can suffer from cumulative effects of sample attrition, and a survey with independent samples, which exhibits a greater sampling error for measuring net changes and trends.
The overlapping panels do however increase the complexity for estimating net changes due to the possibly high correlation. For a statistical production process the estimation can become quite burdensome since estimates are produced not only for the whole survey population but also for many different sub-groupings.
One prominent sample with a rotating panel design is the so called EU Statistics on Income and Living Conditions (EU-SILC), launched 2004. It represents a European standardized survey to generate comparative measures of poverty and social exclusion among the EU Member States and social inclusion policies implemented by Member States are monitored mainly with comparative indicators based on EU-SILC data. The survey is conducted annually in each country with a rotating panel design of 4 waves [1]. In combination with a harmonized survey a set of common indicators, the so called Laeken indicators, were adopted for the countries of the EU [2].
In case of EU-SILC the ultimate cluster approach has been suggested during the second network project for the analysis of EU-SILC (Net-SILC2) as a pragmatic solution for estimating standard errors [3]. This method approximates total sampling variance from the variance between clusters at the first stage, which is called the ultimate cluster approach, see [4]. Standard errors for non-linear indicators are calculated upon linearised variables which have asymptotically the same variance properties as the indicators in question [5]. The approach has been implemented in an R-package called
The survey design and sample size of the EU-SILC is tuned to deliver qualitatively high well-being indicators at national or Nomenclature of Territorial Units for Statistics level 1 (NUTS1), but usually lacks the capability to do the same for smaller subgroups of the survey population, for example estimates on NUTS2 or NUTS3 level.
When it comes to regional indicators based on sample data such as EU-SILC three approaches are commonly known to improve their precision [7]:
Improved size, allocation or design of regional samples; Improved estimation techniques which use auxiliary information; Relaxing requirements of reporting by aggregating information over space, time for indicators.
The first approach requires changes in the data collection which imply considerable cost and needs time for implementation.
Small Area Estimation (SAE) has been developed as an alternative (or complementary) strategy to changes in sample design [8]. Depending on the model specifications and used auxiliary information the methods can yield different results, precision gains and also potential bias. Considering that EU-SILC is a harmonised survey comparing results on an EU-Level can be challenging for users when the underlying methodology can vary greatly. Another model based approach has for example been developed by the World Bank [9], it combines sample data and administrative information and produces regional estimates by mass imputing the variable of interest on a population census.
In this work we present an R package which supports a harmonious approach for estimating standard errors on data collected through a survey with a rotating panel design. The methodology implemented uses bootstrapping techniques for error estimations which, contrary to jack knife replicates, yield consistent estimates for variance even for non-smooth estimators like the median. In order to mimic the panel design of the data requires that sample units must be linked over the panel waves through a unique identifier. Furthermore statistically significant estimates on subgroups of the population can be derived by using multiple consecutive panel waves. From a subject matter perspective this method is however only applicable if the structural patterns for the indicator changes slowly over time.
In a comprehensive study of poverty in NUTS2 regions [10] Statistics Austria has derived standard errors with a bootstrap replication method which partly drew on experience made in the Advanced Methodology for European Laeken Indicators (AMELI) project [11]. Results were compared over several years, indicators and estimation methods. Standard errors of EU-SILC estimates which are cumulated over three successive years were found to be about 25% below that of single years. This implies an increase of effective sample size by approximately 78%. The study in [10] also showed that results from more advanced SAE techniques could generally not provide more stable results.
Methodology
The following subsection presents the used methodology which is applied on multiple consecutive panels of survey data, for instance multiple years of EU-SILC. The methodology contains the following steps, in this order
Draw 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
Rescaled Bootstrap
Bootstrapping has long been around and used widely to estimate confidence intervals and standard errors of point estimates, [12]. Given a random sample
Let
with
sIn context of sample surveys naive bootstrap procedures are not applicable because they do not incorporate the specific sampling design. In the
Single primary sampling units (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. The bootstrap procedure implemented in the Package surveysd automatically detects single PSUs and either merges them together with the strata or cluster with the next least number of PSUs or substitute with the mean value over all bootstrap replicates at the stage which did not contain single PSUs.
Considering panel rotation and split households
For a survey with a rotating penal design the
This procedure is also extended to cases concerning so called split households. A split household occurs when a sampling unit changes address between consecutive panel waves and is followed up on, thus including other cohabitants at the new dress which were not originally in the sample. These new individuals in the split household will inherit the same bootstrap replicate as the sampling unit which changed dress.
Calibration
The uncalibrated bootstrap weights
Variance estimation
Using the calibrated bootstrap sample weights,
with
where
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 one can pool consecutive panels
Section of the UDB data from Austria
Section of the UDB data from Austria
Doing this for all
with
Applying the filter over the time series of estimated
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 [10] 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.
Drawing the bootstrap replicates Calibrating bootstrap replicate weights Estimating standard errors using bootstrap replicate weights
In the following subsections each of these steps is demonstrated using EU-SILC UDB data from Austria and Spain, containing the years 2008 till 2016. The EU-SILC sample for Austria contains a stratified sample design, whereas for Spain a 2-Stage sample design is used. The UDB data does not include all sample design variables and longitudinal identifiers, however they were provided by INE and Statistics Austria to augment the UDB data. Table 1 shows a snippet for the data from Austria, containing identifiers, like hid and pid, regional variables like nuts2 and degree of urbanisation (urban), as well as sociodemographic variables like sex and age groups (ageGroup).
Results from the function calc.stError()
Results from the function calc.stError()
The following sections demonstrate how variance estimates for the amount of people at risk of poverty or social exclusion (AROPE) can be created.
In order to draw bootstrap replicates, information about the sampling design, cluster and strata need to be available. In the case of the UDB data from Austria the stratification variable is equal to NUTS2. To reproduce the rotation for the bootstrap replicates the household id hid need to be fixed over time. First 1000 bootstrap replicates for the UDB data from Austria are drawn
# define stratified 1-Stage cluster sample
datAT_boot
weights
period
In the case of Spain the sampling design is more complex. At the first stage clusters were stratified and sampled from each strata separately. At the the second stage the households were sampled in each previously sampled cluster. For specifying this sampling design the number of clusters, a variable indicating clusters and stratification at the first stage need to be available in the data. Multiple stages can then be supplied using the arguments cluster and strata. The number of sampling units in each stage and strata must also be supplied, in this case total number of clusters (N.cluster) and households (N.households).
# define stratified 2-Stage cluster sample
datES_boot
1000, hid
weights
strata
cluster
totals
period
AROPE standard error for 2015 compared with 3 year average 2014–2016.
Note that the actual number of clusters was not contained in the UDB data but a rough estimate on the number of clusters was derived for the purpose of this demonstration.
With the function
datAT_calib
conP.var
conH.var
Estimating standard errors
Having the calibrated bootstrap replicate weights at hand we are interested in estimating the rate of persons who are at risk of poverty or social exclusion (AROPE), for each year as well as the net change of this estimate from 2016 to 2015 and from 2016 to 2008. The parameters fun and var in
est
fun
period.diff
est
Table 2 shows the results of the calculation. Columns that use the val_ prefix denote the point estimate belonging to the “main weight” of the dataset, which is weight in case of the dataset used here. Columns with the stE_ prefix denote standard errors calculated with bootstrap replicates. The replicates are calculated by using the replicate weights instead of weight when applying the estimator. n denotes the number of observations for the year and N denotes the total weight of those persons.
Grouping
Usually estimates need to be produced not only for each panel wave but also for many different subgroups of the underlying population, for instance for each NUTS2 region, age group, sex and a combination of those. Using the group argument this can easily be supplied to the function calc.stError()
groups
resGroup
“AROPE”, fun
group
period.diff
resGroup
## Calculated point estimates for variable(s)
##
## AROPE
##
## Results hold 297 point estimates for 9 periods in
27 subgroups
##
## Estimated standard error exceeds 10 % of the the
point estimate in 86 cases
From the output we see that we have calculated 297 point estimates for 27 different groupings, which result from
Custom estimators
There are two built-in estimator functions in calc. stError(): weightedSum() and weightedRatio(). In order to define a custom estimator function to be used in fun, the function needs to have two arguments like the example below.
## define custom estimator
myWeightedRatio
}
est2
“AROPE”, fun
period.diff
## [1] TRUE
The parameters x and w can be assumed to be vectors with equal length with w being numeric and x being the column defined in the var argument. It will be called once for each period (in this case year) and for each weight column (in this case weight and the replicate weights).
Comparison with 3-year averages
When producing point estimates for subgroups of the population the underlying sample size might not suffice to produce reliable estimates. As mentioned in the methodological section one can overcome this issue by estimating averages over multiple consecutive panel waves. Using the argument period.mean averages are constructed for each panel wave and grouping, resulting in
groups
resAT
“AROPE”, fun
period.mean
resES
“AROPE”, fun
period.mean
Figure 1 shows the comparison of the standard error of AROPE for each region (nuts2) by degree of urbanisation (urban) when using a single year, 2015, and 3 consecutive years, 2014 to 2016. The dotted line represents the 25% improvement that is to be expected for using the mean over 3 consecutive years [10]. Most estimates are indeed gathering around that line.
Conclusion and outlook
Surveys with a rotating panel design are commonly used when net changes or trends over time of sociodemographic indicators are of interest. The introduction of rotating panels does have a decreasing effect on the variance but results in increasing complexity of the standard error estimates due to the possibly high correlations between panel waves. Using the R package surveysd estimates for standard errors for surveys with a rotating panel design can be derived in a straight forward way. Step by step instructions for the usage of the package are available at
Given the computational flexibility of the replication approach, the
Technical improvement of the replication method may be necessary to represent the actual longitudinal sampling by rotational group. In the present version of the algorithm, bootstrap replicates are generated independently for each wave and only afterwards taken forward for individuals from the same original sample household. Results from bootstrapping and calibration could be displayed in a more userfriendly way. This would require additional classes and possibly some interfaces to the rmarkdown package. Similarly, markup methods for the results of calc.stError() could be added. Currently, there is no way to obtain the actual distribution of bootstrap replicates. Only derived estimates (standard errors and confidence intervals) are returned by calc.stError(). A more flexible approach would be to return the replicates and to define methods that summarize or visualize those replicates.
Footnotes
Acknowledgments
Early parts of the R package surveysd were initially developed in the course of the third Network for the analysis of EU-SILC (Net-SILC3), funded by Eurostat. Authors gratefully acknowledge generous support from several colleagues. Alexander Kowarik, Angelika Meraner and Matthias Till provided fruitfull feedback and extended testing of the softwares.
