Data from complex survey designs require special consideration with regard to estimation of finite population parameters and corresponding variance estimation procedures, as a consequence of significant departures from the simple random sampling assumption. In the past decade a number of statistical software packages have been developed to facilitate the analysis of complex survey data. All these statistical software packages are able to treat samples selected from one sampling frame containing all population units. Dual frame surveys are very useful when it is not possible to guarantee a complete coverage of the target population and may result in considerable cost savings over a single frame design with comparable precision. There are several estimators available in the statistical literature but no existing software covers dual frame estimation procedures. This gap is now filled by package Frames2. In this paper we highlight the main features of the package. The package includes the main estimators in dual frame surveys and also provides interval confidence estimation.
Classic sampling theory usually assumes the existence of one sampling frame containing all finite population units. Then, a probability sample is drawn according to a given sampling design and information collected is used for estimation and inference purposes. In traditional ‘design-based’ inference the population data are regarded as fixed and the randomness comes entirely from the sampling procedure. The most used design-based estimator is the Horvitz-Thompson estimator that is unbiased for the population total if the sampling frame includes all population units, if all sampled units respond and if there is no measurement error. In the presence of auxiliary information, there exist several procedures to obtain more efficient estimators for population means and totals of variables of interest; in particular, customary ratio, regression, raking, post-stratified and calibration estimators. Several software packages have been developed to facilitate the analysis of complex survey data and implement some of these estimators as SAS (SAS Institute Inc. 2013), SPSS (IBM Corporation 2013), Systat (Systat Software Inc. 2009), Stata (Stata Corporation 2015), SUDAAN (Research Triangle Institute 2013) and PCCarp (Fuller et al. 1989). CRAN hosts several R packages that include these design-based methods typically used in survey methodology to treat samples selected from one sampling frame (e.g. survey, Lumley (2014); sampling, Tillé and Matei (2012); laeken, Alfons et al. (2014) or TeachingSampling, Gutierrez Rojas (2014), among others). (Templ 2014) provides a detailed list of packages that includes methods to analyse complex surveys.
In practice, the assumption that the sampling frame contains all population units is rarely met. Often, one finds that sampling from a frame which is known to cover approximately all units in the population is quite expensive while other frames (e.g. special lists of units) are available for cheaper sampling methods. However, the latter usually only cover an unknown or only approximately known fraction of the population. A common example of frame undercoverage is provided by telephone surveys. Estimation could be affected by serious bias due to the lack of a telephone in some households and the generalised use of mobile phones, which are sometimes replacing fixed (land) lines entirely. The potential for coverage error as a result of the exponential growth of the cell-phone only population has led to the development of dual-frame surveys. In these designs, a traditional sample from the landline frame is supplemented with an independent sample from a register of cell-phone numbers.
The dual frame sampling approach assumes that two frames are available
for sampling and that, overall, they cover the entire target population.
The most common situation is the one represented in Figure
1 where the two frames, say frame
Nevertheless, one can face some other situations depending on the
relative positions of the frames. For example, Figure 2
shows the case in which frame
Whatever the scenario, an appropriate choice of the frames results in a better coverage of the target population, which, in turn, leads to a better efficiency of estimators calculated from data from dual frame surveys. This point is particularly important when estimating parameters in rare or elusive populations, where undercoverage errors are usually due to the difficulty of finding individuals showing the characteristic under study when sampling from only one general frame. This issue can be dealt with by incorporating a second frame with a high density of members of the rare population so that the two frames are, together, now complete. Dual frame sampling as a method of improvement of efficiency may seem expensive and unviable, but it is not. In fact, (Hartley 1962) notes that dual frame surveys can result in important cost savings in comparison with single frame ones with a comparable efficiency. As an additional interesting characteristic, dual frame methodology offers the researcher the possibility to consider different data collection procedures and/or different sampling designs, one for each frame. Dual frame surveys have gained much attention and became largely used by statistical agencies and private organizations to take advantage of these benefits.
Standard software packages for complex surveys cannot be used directly when the sample is obtained from a dual frame survey because the classical design-based estimators are severely biased and there is a underestimation of standard errors. Weighted analyses with standard statistical software, with certain modified weights, can yield correct point estimates of population parameters but still yield incorrect results for estimated standard errors. A number of authors have developed methods for estimating population means and totals from dual (or, more generally, multiple) frame surveys but most of these methods require ad-hoc software for their implementation. To the best of our knowledge, there is no software incorporating these estimation procedures for handling dual frame surveys.
Frames2 (Arcos et al. 2015) tries to fill this gap by providing functions for point and interval estimation from dual frame surveys. The paper is organized as follows. In the next section, we provide an overview of the main point estimators proposed so far in the dual frame context and review also jackknife variance estimation as a tool to compare efficiency for all of them. Subsequently, we present package Frames2, discussing guidelines that have been followed to construct it and presenting its principal functions and functionalities. We also provide examples to illustrate how the package works.
Consider again the situation depicted in Figure 1.
Assume we have a finite set of
Let
Let
From data collected in
Different approaches for estimating the population total from dual frame
surveys have been proposed in the literature. (Hartley 1962) suggests
the use of a parameter,
where
is the optimum value for
where hats denote suitable variance and covariance estimators.
(Fuller and Burmeister 1972) introduce information from the estimation of overlap domain size, obtaining the following estimator
where
with
(Bankier 1986) and (Kalton and Anderson 1986) combine all sampling units coming from
the two frames,
or, summarizing,
Hence, the estimator can be expressed in the form
with
where
When frame sizes,
where
(Skinner and Rao 1996) use a pseudo maximum likelihood approach to extend to complex designs the maximum likelihood estimator proposed by (Fuller and Burmeister 1972) only for simple random sampling without replacement. The resulting estimator is given by
where
minimizes the variance of
where, in this case,
More recently, (Rao and Wu 2010) proposed a pseudo empirical likelihood estimator for the population mean based on poststratified samples. The estimator is computed as
where, in this case,
Recently, (Ranalli et al. 2013) extended calibration procedures to estimation
from dual frame sampling assuming that some kind of auxiliary
information is available. For example, assuming there are
where weights
where
Then, with a similar approach to that of
with weights
An estimator of the variance of any calibration estimator can be obtained using Deville’s method (Deville 1993) through the following expression
where
Some of the estimators described above are particular types of
calibration estimators. For example, estimator ((9)) can be
obtained as a particular case of "raking"
method is
selected for calibration. Having noted this, one can use
((17)) to calculate an estimator of variance of
((9)). See (Ranalli et al. 2013) for more details.
Table 1 shows a summary of the previous dual frame
estimators according to the auxiliary information required. It can be
noted that Hartley, FB and BKA estimators can be computed even when no
information is available, but they cannot incorporate some auxiliary
information when available. PML and SFRR can incorporate information on
None | known | |||||||
Hartley | ||||||||
FB | ||||||||
PML | ||||||||
PEL | ||||||||
CalDF | ||||||||
BKA* | ||||||||
SFRR* | ||||||||
CalSF* |
(*) Inclusion probabilities are known in overlap domain
Variance estimation methods exposed so far depend on each specific estimator, so comparisons between variance estimations may lead to incorrect conclusions. Instead, one can consider jackknife, originally proposed by (Quenouille 1949, 1956) (see Wolter (2007) for a detailed description of this method in survey sampling) and extended to dual frame surveys by (Lohr and Rao 2000), which can be used to estimate variances irrespective of the type of estimator allowing us to compare estimated efficiency for different estimators.
For a non-stratified design in each frame, the jackknife estimator of
the variance for any of the estimators described, generically denoted by
with
Jackknife may present an important bias when designs are without
replacement. One could, then, incorporate an approximate
finite-population correction to estimation to achieve unbiasedness. For
example, assuming that a finite-population correction is needed in frame
Consider now a stratified design in each frame, where frame
where
Stratified cluster sampling is very common in practice. Now we
illustrate the jackknife estimator when a stratified sample of clusters
is selected. Suppose that frame
To define the jackknife estimator of the variance, let
where
Frames2 is a new R package for point and interval estimation from dual
frame sampling. It consists of eight main functions (Hartley
, FB
,
BKA
, SFRR
, PML
, PEL
, CalSF
and CalDF
), each of them
implementing one of the estimators described in the previous sections.
The package also includes an additional function called Compare
which
provides a summary with all possible estimators that can be computed
from the information provided as input. Moreover, six extra functions
implementing auxiliary operations, like computation of Horvitz-Thompson
estimators or of the covariance between two Horvitz-Thompson estimators,
have also been included in the package to achieve a more understandable
code. Finally, the package includes eight more functions, one for each
estimator, for the calculation of confidence intervals based on the
jackknife variance estimator.
A notable feature of these functions is the strong argument check.
Functions check general aspects such as the presence of NA
or NaN
values in its arguments, the number of main variables considered in the
frames (that should match), the length of the arguments in each frame
(that should also match) and the values for arguments indicating the
domain each unit belongs to (which only can be "a"
or "ab"
for frame
"b"
or "ba"
for frame
Much attention has also been devoted to computational efficiency. Frequently, populations in a survey are extremely large or it is needed to keep sampling error below a certain value. As a consequence, one needs to consider large sample sizes, often in the order of tens of thousands sampling units. In these situations, computational efficiency of functions is essential, particularly when several variables are considered. Otherwise, the user can face high runtimes and heavy computational loads. In this sense, functions of Frames2 are developed according to strict efficiency measures, using the power of R to use matrix calculation to avoid loops and increase the computational efficiency.
Table 2 shows user and system times necessary to compute estimators using an Intel(R) Core(TM) i7-3770 at 3.40 GHz when different sample sizes are considered. Elapsed time is also included to get an idea about the real time user needs to get estimations.
user | system | elapsed | ||||
Hartley | 0.01 | 0.02 | 0.04 | |||
FB | 0.05 | <0.01 | 0.07 | |||
BKA | 0.03 | <0.01 | 0.05 | |||
PML | 0.02 | 0.02 | 0.03 | |||
SFRR | 0.03 | 0.03 | 0.07 | |||
CalSF | 0.03 | <0.01 | 0.06 | |||
CalDF | 0.04 | 0.01 | 0.05 | |||
Hartley | 0.11 | 0.06 | 0.19 | |||
FB | 0.27 | 0.07 | 0.32 | |||
BKA | 0.13 | 0.05 | 0.17 | |||
PML | 0.16 | 0.02 | 0.18 | |||
SFRR | 0.42 | 0.12 | 0.54 | |||
CalSF | 0.20 | 0.08 | 0.30 | |||
CalDF | 0.22 | 0.07 | 0.31 |
Functions of Frames2 have been implemented from an user-oriented perspective to increase usability. In this sense, most input parameters (which are the communication channel between the user and the function) are divided into two groups, depending on the frame they come from. This is to adapt functions as much as possible to the usual estimation procedure, in which the first step is to draw two independent samples, one from each frame. On the other hand, estimation details are managed internally by functions so that they are not visible for the user, who does not need to manage them.
Construction of functions has been carried out so that they perform properly in as many situations as possible. As noted in the introductory section, one can face several situations when using two sampling frames depending on their relative positions. Although the most common is the one depicted in Figure 1, cases shown in Figures 2 and 3 may arise as well. All estimators described except PEL can be modified to cover these three situations, so corresponding functions of Frames2 include necessary changes to produce estimates irrespective of the situation.
On the other hand, it is usual, when conducting a survey, to collect
information on many variables of interest. To adapt to such situations,
all functions are programmed to produce estimates when there are more
than one variable of interest with only one call. To this end,
parameters containing information about main variables observed in each
frame can be either vectors, when only one variable is considered or
matrices or data frames, when there are several variables under study.
Cases in which the main aim of the survey is the estimation of
population means or proportions are also very frequent. Hence, from the
estimation of the population total for a variable, functions compute
estimation of the mean as
To get maximum flexibility, functions have been programmed to calculate
estimates in cases in which the user disposes of first and second order
inclusion probabilities and in those other in which only first order
ones are available, indistinctly. Knowledge of both first and second
order inclusion probabilities is a strong assumption that does not
always occur in practice. However, when calculating most of the
estimators described in previous sections, second order inclusion
probabilities are needed in many steps of the estimation procedure,
mainly in computing estimated variances of a Horvitz-Thompson estimator
or estimated covariances between two Horvitz-Thompson estimators. As an
alternative, one can obtain variance estimations from only first order
inclusion probabilities applying Deville’s method reported in
((17)), by substituting residuals
As can be deduced from previous sections, an essential aspect when
computing estimates in dual frames is to know the domain each unit
belongs to. Character vectors domains_A
and domains_B
are used for
this purpose. The former can take values "a"
or "ab"
, while the
latter can take values "b"
or "ba"
. Any other value will be
considered as incorrect.
To illustrate how functions operate, we use the data sets DatA
and
DatB
, both included in the package. DatA
contains information about
Both data sets contain information about the same variables. To better
understand their structure, we report the first three rows of DatA
:
> library(Frames2)
> data(DatA)
> head(DatA, 3)
Domain Feed Clo Lei Inc Tax M2 Size ProbA ProbB Stratum
1 a 194.48 38.79 23.66 2452.07 112.90 0.00 0 0.02063274 0.0000000 1
2 a 250.23 16.92 22.68 2052.37 106.99 0.00 0 0.02063274 0.0000000 1
3 ab 199.95 24.50 23.24 2138.24 121.16 127.41 2 0.02063274 0.1133501 1
Each data set incorporates information about three main variables:
Feeding, Clothing and Leisure. Additionally, there are two auxiliary
variables for the units in frame
Numerical square matrices PiklA
and PiklB
, with dimensions
PiklA
is shown.
> data(PiklA)
> PiklA[1:6, 1:6]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.020632737 0.000397876 0.000397876 0.000397876 0.000397876 0.000397876
[2,] 0.000397876 0.020632737 0.000397876 0.000397876 0.000397876 0.000397876
[3,] 0.000397876 0.000397876 0.020632737 0.000397876 0.000397876 0.000397876
[4,] 0.000397876 0.000397876 0.000397876 0.020632737 0.000397876 0.000397876
[5,] 0.000397876 0.000397876 0.000397876 0.000397876 0.020632737 0.000397876
[6,] 0.000397876 0.000397876 0.000397876 0.000397876 0.000397876 0.020632737
When there is no further information than the one on the variables of interest, one can calculate some of the estimators described in previous sections (as, for example, ((1)) or ((4))) as follows
> library(Frames2)
>
> data(DatA)
> data(DatB)
> data(PiklA)
> data(PiklB)
>
> yA <- with(DatA, data.frame(Feed, Clo))
> yB <- with(DatB, data.frame(Feed, Clo))
>
> ## Estimation for variables Feeding and Clothing using Hartley and
> ## Fuller-Burmeister estimators with first and second order probabilities known
> Hartley(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain)
Estimation:
Feed Clo
Total 586959.9820 71967.62214
Mean 246.0429 30.16751
> FB(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain)
Estimation:
Feed Clo
Total 591665.5078 72064.99223
Mean 248.0153 30.20832
> ## This is how estimates change when only first order probabilities are considered
> Hartley(yA, yB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain)
Estimation:
Feed Clo
Total 570867.8042 69473.86532
Mean 247.9484 30.17499
> FB(yA, yB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain)
Estimation:
Feed Clo
Total 571971.9511 69500.11448
Mean 248.4279 30.18639
As the result, an object of class "EstimatorDF"
is returned, showing,
by default, estimations for the population total and mean for the two
considered variables. In general, summary
.
> summary(Hartley(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain))
Call:
Hartley(ysA = yA, ysB = yB, pi_A = PiklA, pi_B = PiklB, domains_A = DatA$Domain,
domains_B = DatB$Domain)
Estimation:
Feed Clo
Total 586959.9820 71967.62214
Mean 246.0429 30.16751
Variance Estimation:
Feed Clo
Var. Total 2.437952e+08 4.728875e+06
Var. Mean 4.283804e+01 8.309261e-01
Total Domain Estimations:
Feed Clo
Total dom. a 263233.1 31476.84
Total dom. ab 166651.7 21494.96
Total dom. b 164559.2 20451.85
Total dom. ba 128704.7 15547.49
Mean Domain Estimations:
Feed Clo
Mean dom. a 251.8133 30.11129
Mean dom. ab 241.6468 31.16792
Mean dom. b 242.2443 30.10675
Mean dom. ba 251.5291 30.38466
Parameters:
Feed Clo
theta 0.8027766 0.7551851
The previous output shows in the component Estimation
the estimations
of the population total and the population mean computed using the
Harley estimator, that is, Variance Estimation
. In the section Total Domain Estimations
we can
see estimations Mean Domain Estimations
. Finally,
This additional information depends on the way each estimator is
formulated. Thus, for example, extra information will include a
parameter component when applied to a call to the Fuller-Burmeister
estimator (and values of estimates for
Results slightly change when a confidence interval is required. In that
case, the user has to indicate the confidence level desired for the
interval through argument conf_level
(default is NULL
) and add it to
the list of input parameters. The function calculates, then, a
confidence interval based on the pivotal method. This method yields a
confidence interval as follows:
> Hartley(yA, yB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain, 0.95)
Estimation and 95 % Confidence Intervals:
Feed Clo
Total 570867.8042 69473.86532
Lower Bound 511904.6588 61756.37677
Upper Bound 629830.9496 77191.35387
Mean 247.9484 30.17499
Lower Bound 222.3386 26.82301
Upper Bound 273.5582 33.52697
> FB(yA, yB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain, 0.95)
Estimation and 95 % Confidence Intervals:
Feed Clo
Total 571971.9511 69500.11448
Lower Bound 513045.7170 61802.57411
Upper Bound 630898.1852 77197.65484
Mean 248.4279 30.18639
Lower Bound 222.8342 26.84307
Upper Bound 274.0217 33.52971
For estimators constructed as ((7)), numeric vectors
pik_ab_B
and pik_ba_A
of lengths pik_ab_B
represents first order inclusion
probabilities according to sampling design in frame pik_ba_A
contains first order inclusion probabilities according
to the sampling design in frame
> yA <- with(DatA, data.frame(Feed, Clo, Lei))
> yB <- with(DatB, data.frame(Feed, Clo, Lei))
>
> ## Bankier-Kalton-Anderson estimation and a 95% confidence
> ## interval for the three main variables
> BKA(yA, yB, DatA$ProbA, DatB$ProbB, DatA$ProbB, DatB$ProbA, DatA$Domain,
+ DatB$Domain, 0.95)
Estimation and 95 % Confidence Intervals:
Feed Clo Lei
Total 566434.3200 68959.26705 50953.07583
Lower Bound 624569.2139 76538.11015 56036.23578
Upper Bound 508299.4262 61380.42395 45869.91588
Mean 247.8845 30.17814 22.29822
Lower Bound 273.3257 33.49482 24.52273
Upper Bound 222.4434 26.86147 20.07372
Note that these examples include just a few of the estimators that can
be used when no auxiliary information is known. As noted in Table
1, other estimators, as those in ((13)),
((14)) or ((16)), can be also calculated in this
case. In this context, function Compare
is quite useful, since it
returns all possible estimators that can be computed according to the
information provided as input.
> Compare(yA, yB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain)
$Hartley
Estimation:
Feed Clo Lei
Total 570867.8042 69473.86532 51284.2727
Mean 247.9484 30.17499 22.2746
$FullerBurmeister
Estimation:
Feed Clo Lei
Total 571971.9511 69500.11448 51210.03819
Mean 248.4279 30.18639 22.24236
$PEL
Estimation:
Feed Clo Lei
Total 1.791588e+08 2.663164e+06 1.455533e+06
Mean 2.479314e+02 3.011373e+01 2.235969e+01
$Calibration_DF
Estimation:
Feed Clo Lei
Total 595162.2604 72214.13351 53108.5059
Mean 248.8422 30.19332 22.2051
Using appropriate indicator variables as variables of interest, one can also estimate the overlap domain size, as shown below:
> indA <- as.integer(DatA$Domain == "ab")
> indB <- as.integer(DatB$Domain == "ba")
>
> Hartley(indA, indB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain)
Estimation:
[,1]
Total 534.2743208
Mean 0.2320545
> BKA(indA, indB, DatA$ProbA, DatB$ProbB, DatA$ProbB, DatB$ProbA, DatA$Domain,
+ DatB$Domain)
Estimation:
[,1]
Total 560.4121771
Mean 0.2452491
For estimators requiring frame sizes known, as ((9)) or
((10)), it is needed to incorporate two additional input
arguments, N_A
and N_B
. There is also a group of estimators,
including ((13)) and ((16)), that, even though able
to provide estimations without the need of auxiliary information, can
use frame sizes to improve their precision. The following examples show
the performance of these estimators.
> ## SFRR estimator and CalSF estimator with frame sizes as auxiliary information
> ## using method "raking" for the calibration for the three main variables
> SFRR (yA, yB, DatA$ProbA, DatB$ProbB, DatA$ProbB, DatB$ProbA, DatA$Domain,
+ DatB$Domain, N_A = 1735, N_B = 1191)
Estimation:
Feed Clo Lei
Total 584713.4070 71086.18669 52423.74035
Mean 248.2219 30.17743 22.25487
> CalSF(yA, yB, DatA$ProbA, DatB$ProbB, DatA$ProbB, DatB$ProbA, DatA$Domain,
+ DatB$Domain, N_A = 1735, N_B = 1191, met = "raking")
Estimation:
Feed Clo Lei
Total 584713.4070 71086.18669 52423.74035
Mean 248.2219 30.17743 22.25487
As highlighted previously, both results match. Note that the argument
met
of the SF calibration estimator indicates the method used in the
calibration procedure. The possibility of choosing the calibration
method is given by the fact that computation of both SF and DF
calibration estimators is based on the function calib
from package
sampling (Tillé and Matei 2012), which can manage three different calibration
methods, each one associated with one particular distance measure. These
methods are: linear, raking and logit.
The requirement of knowing probabilities of inclusion in both frames for the units in the overlap domain may be restrictive in some cases. As an alternative, in cases where frame sizes are known but this condition is not met, it is possible to calculate dual frame estimators as ((10)), ((13)) or ((14)). Next, it is illustrated how to obtain some of these estimators with Frames2.
> ## Estimates for the three main variables using PML, PEL and CalDF
> ## with frame sizes as auxiliary information in PEL and CalDF
> PML(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191)
Estimation:
Feed Clo Lei
Total 593085.4467 72272.73759 53287.68044
Mean 248.0966 30.23277 22.29104
> PEL(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191)
Estimation:
Feed Clo Lei
Total 590425.4843 72211.61334 53258.38286
Mean 247.4958 30.26982 22.32497
> CalDF(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191)
Estimation:
Feed Clo Lei
Total 587502.4374 71368.45308 52490.98852
Mean 248.7193 30.21385 22.22207
To calculate the PEL
estimator, computational algorithms for the
pseudo empirical likelihood method for the analysis of complex survey
data presented by (Wu 2005) are used.
In addition to the frame sizes, in some cases, it is possible to know
the size of the overlap domain, N_ab
, as shown
below
> ## Estimates for the three main variables using PEL estimator
> ## with frame sizes and overlap domain size as auxiliary information
> PEL(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191,
+ N_ab = 601)
Estimation:
Feed Clo Lei
Total 575289.2186 70429.95642 51894.32490
Mean 247.4362 30.29245 22.32014
> ## Calibration estimators with the same auxiliary information
> ## Estimates do not change when raking method is used for the calibration
> CalSF(yA, yB, PiklA, PiklB, DatA$ProbB, DatB$ProbA, DatA$Domain, DatB$Domain,
+ N_A = 1735, N_B = 1191, N_ab = 601)
Estimation:
Feed Clo Lei
Total 577163.6066 70173.20412 51726.19862
Mean 248.2424 30.18202 22.24783
> CalSF(yA, yB, PiklA, PiklB, DatA$ProbB, DatB$ProbA, DatA$Domain, DatB$Domain,
+ N_A = 1735, N_B = 1191, N_ab = 601, met = "raking")
Estimation:
Feed Clo Lei
Total 577163.6067 70173.20414 51726.19863
Mean 248.2424 30.18202 22.24783
> CalDF(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191,
+ N_ab = 601)
Estimation:
Feed Clo Lei
Total 578691.1756 70246.32319 51600.78973
Mean 248.8994 30.21347 22.19389
> CalDF(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191,
+ N_ab = 601, met = "raking")
Estimation:
Feed Clo Lei
Total 578691.1763 70246.32328 51600.78979
Mean 248.8994 30.21347 22.19389
Note that in this case, calibration estimators provide the same results irrespective of the distance function employed. This is an interesting property that calibration estimators show only in the case in which all the domain sizes are known and used for calibration (see Deville (1993)).
On the other hand, some of the estimators are defined such that they can
incorporate auxiliary information into the estimation process. This is
the case for estimators ((13)), ((14)) and
((16)). Functions implementing them are also able to manage
auxiliary information. To achieve maximum flexibility, functions
implementing estimators ((13)), ((14)) and
((16)) are prepared to deal with auxiliary information when
it is available only in frame xsAFrameA
and
xsBFrameA
, numeric vectors, matrices or data frames (depending on the
number of auxiliary variables in the frame); and XA
, a numeric value
or vector of length indicating population totals for the auxiliary
variables considered in frame xsAFrameB
, xsBFrameB
and XB
. If auxiliary information is available
in the whole population, this must be indicated through parameters xsT
and X
. In the following example, one can see how to calculate
estimators using different type of auxiliary information
> ## PEL, CalSF and CalDF estimators for the three main variables
> ## using Income as auxiliary variable in frame A and Metres2 as auxiliary
> ## variable in frame B assuming frame sizes known
> PEL(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191,
+ xsAFrameA = DatA$Inc, xsBFrameA = DatB$Inc, xsAFrameB = DatA$M2,
+ xsBFrameB = DatB$M2, XA = 4300260, XB = 176553)
Estimation:
Feed Clo Lei
Total 587742.7193 71809.56826 53094.20112
Mean 246.3713 30.10129 22.25614
> CalSF(yA, yB, PiklA, PiklB, DatA$ProbB, DatB$ProbA, DatA$Domain, DatB$Domain,
+ N_A = 1735, N_B = 1191, xsAFrameA = DatA$Inc, xsBFrameA = DatB$Inc,
+ xsAFrameB = DatA$M2, xsBFrameB = DatB$M2, XA = 4300260, XB = 176553)
Estimation:
Feed Clo Lei
Total 582398.3181 70897.88438 52252.24741
Mean 247.5819 30.13922 22.21282
> CalDF(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191,
+ xsAFrameA = DatA$Inc, xsBFrameA = DatB$Inc, xsAFrameB = DatA$M2,
+ xsBFrameB = DatB$M2, XA = 4300260, XB = 176553)
Estimation:
Feed Clo Lei
Total 585185.4497 71194.61148 52346.43878
Mean 247.8075 30.14866 22.16705
> ## Now, assume that overlap domain size is also known
> PEL(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191,
+ N_ab = 601, xsAFrameA = DatA$Inc, xsBFrameA = DatB$Inc,
+ xsAFrameB = DatA$M2, xsBFrameB = DatB$M2, XA = 4300260, XB = 176553)
Estimation:
Feed Clo Lei
Total 572611.6997 69991.74803 51737.56089
Mean 246.2846 30.10398 22.25271
> CalSF(yA, yB, PiklA, PiklB, DatA$ProbB, DatB$ProbA, DatA$Domain, DatB$Domain,
+ N_A = 1735, N_B = 1191, N_ab = 601, xsAFrameA = DatA$Inc, xsBFrameA = DatB$Inc,
+ xsAFrameB = DatA$M2, xsBFrameB = DatB$M2, XA = 4300260, XB = 176553)
Estimation:
Feed Clo Lei
Total 575636.7876 70076.78485 51628.27583
Mean 247.5857 30.14055 22.20571
> CalDF(yA, yB, PiklA, PiklB, DatA$Domain, DatB$Domain, N_A = 1735, N_B = 1191,
+ N_ab = 601, xsAFrameA = DatA$Inc, xsBFrameA = DatB$Inc,
+ xsAFrameB = DatA$M2, xsBFrameB = DatB$M2, XA = 4300260, XB = 176553)
Estimation:
Feed Clo Lei
Total 576630.7609 70102.0037 51477.16737
Mean 248.0132 30.1514 22.14072
Finally, eight additional functions have been included, each of them
calculating confidence intervals based on jackknife variance estimation
for each estimator. To carry out variance estimation using the jackknife
method, in addition to parameters to calculate each specific estimator,
the user has to indicate through arguments sdA
and sdB
the sampling
design applied in each frame. Possible values are "srs"
(simple random
sampling without replacement), "str"
(stratified sampling), "pps"
(probabilities proportional to size sampling), "clu"
(cluster
sampling) or "strclu"
(stratified cluster sampling). Default is
"srs"
for both frames. If a stratified or a cluster sampling has been
carried out in one of the frames, it is needed to include information
about the strata or the clusters. Furthermore, the user is able to
include a finite population correction factor in each frame by setting
to TRUE
the parameters fcpA
and fcpB
, set by default to FALSE
.
As the main purpose of the functions is to obtain confidence intervals,
parameter conf_level
is now mandatory. As noted, these functions can
be used, for example, to make comparisons between efficiency of
estimators, as shown in the next example.
> ## Confidence intervals through jackknife for the three main variables
> ## for estimators defined under the so called single frame approach with
> ## a stratified random sampling in frame A and a simple random sampling
> ## without replacement in frame B. Finite population correction factor
> ## is required for frame A
> JackBKA (yA, yB, DatA$ProbA, DatB$ProbB, DatA$ProbB, DatB$ProbA, DatA$Domain,
+ DatB$Domain, conf_level = 0.95, sdA = "str", strA = DatA$Stratum, fcpA = TRUE)
Feed Clo Lei
Total 566434.3200 68959.26705 50953.07583
Jack Upper End 610992.1346 74715.89841 54717.32664
Jack Lower End 521876.5055 63202.63570 47188.82502
Mean 247.8845 30.17814 22.29822
Jack Upper End 267.3840 32.69738 23.94555
Jack Lower End 228.3850 27.65891 20.65090
> JackSFRR(yA, yB, DatA$ProbA, DatB$ProbB, DatA$ProbB, DatB$ProbA, DatA$Domain,
+ DatB$Domain, N_A = 1735, N_B = 1191, conf_level = 0.95, sdA = "str",
+ strA = DatA$Stratum, fcpA = TRUE)
Feed Clo Lei
Total 584713.4070 71086.18669 52423.74035
Jack Upper End 619959.0338 76576.74587 55204.67760
Jack Lower End 549467.7802 65595.62751 49642.80309
Mean 248.2219 30.17743 22.25487
Jack Upper End 263.1843 32.50828 23.43543
Jack Lower End 233.2595 27.84659 21.07431
> JackCalSF(yA, yB, DatA$ProbA, DatB$ProbB, DatA$ProbB, DatB$ProbA, DatA$Domain,
+ DatB$Domain, N_A = 1735, N_B = 1191, N_ab = 601, conf_level = 0.95, sdA = "str",
+ strA = DatA$Stratum, fcpA = TRUE)
Feed Clo Lei
Total 577163.6066 70173.20412 51726.19862
Jack Upper End 599105.4275 73516.53187 53165.97439
Jack Lower End 555221.7858 66829.87636 50286.42285
Mean 248.2424 30.18202 22.24783
Jack Upper End 257.6798 31.62001 22.86709
Jack Lower End 238.8051 28.74403 21.62857
> ## Same for a selection of dual frame estimators
> JackHartley (yA, yB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain,
+ conf_level = 0.95, sdA = "str", strA = DatA$Stratum, fcpA = TRUE)
Feed Clo Lei
Total 570867.8042 69473.86532 51284.27265
Jack Upper End 610664.7131 74907.33129 54782.33083
Jack Lower End 531070.8954 64040.39934 47786.21447
Mean 247.9484 30.17499 22.27460
Jack Upper End 265.2336 32.53494 23.79393
Jack Lower End 230.6631 27.81504 20.75527
> JackPML(yA, yB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain,
+ N_A = 1735, N_B = 1191, conf_level = 0.95, sdA = "str", strA = DatA$Stratum,
+ fcpA = TRUE)
Feed Clo Lei
Total 594400.6320 72430.05834 53408.30337
Jack Upper End 626443.7529 76885.06491 56003.77592
Jack Lower End 562357.5111 67975.05176 50812.83082
Mean 248.0934 30.23115 22.29178
Jack Upper End 261.4677 32.09060 23.37509
Jack Lower End 234.7191 28.37171 21.20847
> JackCalDF(yA, yB, DatA$ProbA, DatB$ProbB, DatA$Domain, DatB$Domain, N_A = 1735,
+ N_B = 1191, N_ab = 601, conf_level = 0.95, sdA = "str", strA = DatA$Stratum,
+ fcpA = TRUE)
Feed Clo Lei
Total 578895.6961 70230.11306 51570.55683
Jack Upper End 601626.7000 73614.66702 53037.42260
Jack Lower End 556164.6921 66845.55910 50103.69107
Mean 248.9874 30.20650 22.18088
Jack Upper End 258.7642 31.66222 22.81179
Jack Lower End 239.2106 28.75078 21.54997
In the example above data are separated into two data sets DatA
and
DatB
containing domain information. But in practice, it is common to
have a joint data set including units from both samples in which there
is not a specific variable indicating the domain where each individual
is placed. However, we can easily split the dataset and format it, so
functions of Frames2 can be applied. To illustrate how to do this, we
are going to use dataset Dat
, which includes some of the variables
collected in a real dual frame survey.
Data included in Dat
comes from an opinion survey on the Andalusian
population with respect to immigration. This survey is conducted using
telephone interviews on adults using two sampling frames: one for
landlines and another one for cell phones. From the landline frame, a
stratified sample of size 1919 was drawn, while from the cell phone
frame, a sample of size 483 is drawn using simple random sampling
without replacement. First-order inclusion probabilities were computed
from a stratified random design in the landline frame and modified
taking into account the number of fixed lines and adults in the
household. In the cell phone frame first-order inclusion probabilities
were computed and modified, given the number of cell phone numbers per
individual. At the time of data collection, frame sizes of land and cell
phones were 4,982,920 and 5,707,655, respectively, and the total
population size was 6,350,916.
The data set includes information about 7 variables: Drawnby
, which
takes value 1 if the unit comes from the landline sample and value 2 if
it comes from the cell phone sample; Stratum
, which indicates the
stratum each unit belongs to (for individuals in the cell phone frame,
value of this variable is NA
); Opinion
the response to the question:
“Do you think that immigrants currently living in Andalusia are quite a
lot?” with value 1 representing “yes” and value 0 representing “no”;
Landline
and Cell
, which record whether the unit possess a landline
or a cell phone, respectively. First order inclusion probabilities are
also included in the data set.
> data(Dat)
> head(Dat, 3)
Drawnby Stratum Opinion Landline Cell ProbLandline ProbCell
1 1 2 0 1 1 0.000673623 8.49e-05
2 1 5 1 1 1 0.002193297 5.86e-05
3 1 1 0 1 1 0.001831489 7.81e-05
From the data of this survey we wish to estimate the number of people in
Andalusia thinking that immigrants currently living in this region are
quite a lot. In order to use functions of package Frames2, we need to
split this dataset. The variables we will use to do this are Drawnby
and Landline
and Cell
.
> attach(Dat)
> ## We can split the original dataset in four new different
> ## datasets, each one corresponding to one domain.
>
> DomainOnlyLandline <- Dat[Landline == 1 & Cell == 0,]
> DomainBothLandline <- Dat[Drawnby == 1 & Landline == 1 & Cell == 1,]
> DomainOnlyCell <- Dat[Landline == 0 & Cell == 1,]
> DomainBothCell <- Dat[Drawnby == 2 & Landline == 1 & Cell == 1,]
>
> ## From the domain datasets, we can build frame datasets
>
> FrameLandline <- rbind(DomainOnlyLandline, DomainBothLandline)
> FrameCell <- rbind(DomainOnlyCell, DomainBothCell)
>
> ## Finally, we only need to label domain of each unit using "a", "b",
> ## "ab" or "ba"
>
> Domain <- c(rep("a", nrow(DomainOnlyLandline)), rep("ab", nrow(DomainBothLandline)))
> FrameLandline <- cbind(FrameLandline, Domain)
>
> Domain <- c(rep("b", nrow(DomainOnlyCell)), rep("ba", nrow(DomainBothCell)))
> FrameCell <- cbind(FrameCell, Domain)
Now dual frame estimators, as PML estimator, can be computed:
> summary(PML(FrameLandline$Opinion, FrameCell$Opinion, FrameLandline$ProbLandline,
+ FrameCell$ProbCell, FrameLandline$Domain, FrameCell$Domain, N_A = 4982920,
+ N_B = 5707655))
Call:
PML(ysA = FrameLandline$Opinion, ysB = FrameCell$Opinion,
pi_A = FrameLandline$ProbLandline, pi_B = FrameCell$ProbCell,
domains_A = FrameLandline$Domain, domains_B = FrameCell$Domain,
N_A = 4982920, N_B = 5707655)
Estimation:
[,1]
Total 3.231325e+06
Mean 4.635634e-01
Variance Estimation:
[,1]
Var. Total 1.784362e+10
Var. Mean 3.672317e-04
Total Domain Estimations:
[,1]
Total dom. a 219145.1
Total dom. ab 2318841.9
Total dom. b 1346646.1
Total dom. ba 1457501.0
Mean Domain Estimations:
[,1]
Mean dom. a 0.4438149
Mean dom. ab 0.4990548
Mean dom. b 0.4172797
Mean dom. ba 0.4674919
Parameters:
gamma 0.3211534
As the overlap domain size is known, we can include additionally this information in the process and compute more accurate estimators as CalDF and CalSF.
> summary(CalDF(FrameLandline$Opinion, FrameCell$Opinion, FrameLandline$ProbLandline,
+ FrameCell$ProbCell, FrameLandline$Domain, FrameCell$Domain, N_A = 4982920,
+ N_B = 5707655, N_ab = 4339659))
Call:
CalDF(ysA = FrameLandline$Opinion, ysB = FrameCell$Opinion,
pi_A = FrameLandline$ProbLandline, pi_B = FrameCell$ProbCell,
domains_A = FrameLandline$Domain, domains_B = FrameCell$Domain,
N_A = 4982920, N_B = 5707655, N_ab = 4339659)
Estimation:
[,1]
Total 2.985028e+06
Mean 4.700153e-01
Variance Estimation:
[,1]
Var. Total 1.478990e+10
Var. Mean 3.666844e-04
Parameters:
eta 0.7296841
>
> summary(CalSF(FrameLandline$Opinion, FrameCell$Opinion, FrameLandline$ProbLandline,
+ FrameCell$ProbCell, FrameLandline$ProbCell, FrameCell$ProbLandline,
+ FrameLandline$Domain, FrameCell$Domain, N_A = 4982920, N_B = 5707655,
+ N_ab = 4339659))
Call:
CalSF(ysA = FrameLandline$Opinion, ysB = FrameCell$Opinion,
pi_A = FrameLandline$ProbLandline, pi_B = FrameCell$ProbCell,
pik_ab_B = FrameLandline$ProbCell, pik_ba_A = FrameCell$ProbLandline,
domains_A = FrameLandline$Domain, domains_B = FrameCell$Domain, N_A = 4982920,
N_B = 5707655, N_ab = 4339659)
Estimation:
[,1]
Total 2.986787e+06
Mean 4.702923e-01
Variance Estimation:
[,1]
Var. Total 1.442969e+10
Var. Mean 3.577539e-04
Observe that the greater the information included in the estimation process is, the greater is the accuracy of the estimates.
The statistical literature about dual frame surveys started around 1960 and its development has evolved very quickly because these surveys are largely used by statistical agencies and private organizations to decrease sampling costs and to reduce frame undercoverage errors that could occur with the use of a single sampling frame.
Dual frame surveys can be more complicated to design and more complicated to analyze than those that use one frame only. There are several estimators of the population total available in the statistical literature. These estimators rely on weight adjustments to compensate for the multiplicity of the units in the overlap domain. Some of these estimators allow the handling of different types of auxiliary information at different levels. Nevertheless, none of the existing statistical software implements all of these estimators.
In this article we illustrate Frames2, a new R package for point and interval estimation in dual frame context. Functions in the package implement the most important estimators in the literature for population totals and means. We include two procedures (Pseudo-Empirical-Likelihood approach and calibration approach) to incorporate auxiliary information about frame sizes and also about one or several auxiliary variables in one or two frames. Post-stratification, raking ratio or regression estimation are all encompassed as particular cases of these estimation procedures. Additional functions for confidence interval estimation based on the jackknife variance estimation have been included as well.
The functionalities of the package Frames2 have been illustrated using
several data sets DatA
, DatB
and Dat
(included in the package)
corresponding to different complex surveys. We envision future additions
to the package that will allow for extensions to more than two frames.
Finally, we would like to direct the reader to the package vignettes
named "estimation"
(Estimation in a dual frame context) and
"formatting.data"
(Splitting and formatting data in a dual frame
context) for further examples and background information.
This study was partially supported by Ministerio de Educación, Cultura y Deporte (grant MTM2012-35650 and FPU grant program, Spain), by Consejerı́a de Economı́a, Innovación, Ciencia y Empleo (grant SEJ2954, Junta de Andalucı́a, Spain) and under the support of the project PRIN-SURWEY (grant 2012F42NS8, Italy). The authors are grateful to Manuel Trujillo, (IESA, Institute of Advanced Social Studies) for providing data and information about the OPIA (Opinions and Attitudes of the Andalusian Population regarding Immigration) survey. This paper and software have been considerably improved by the comments of the editor and the reviewers.
survey, sampling, laeken, TeachingSampling, Frames2
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Arcos, et al., "Frames2: A Package for Estimation in Dual Frame Surveys", The R Journal, 2015
BibTeX citation
@article{RJ-2015-005, author = {Arcos, Antonio and Molina, David and Rannalli, Maria Giovanna and Rueda, María del Mar}, title = {Frames2: A Package for Estimation in Dual Frame Surveys}, journal = {The R Journal}, year = {2015}, note = {https://rjournal.github.io/}, volume = {7}, issue = {1}, issn = {2073-4859}, pages = {52-72} }