PracTools : Computations for Design of Finite Population Samples

PracTools is an R package with functions that compute sample sizes for various types of finite population sampling designs when totals or means are estimated. One-, two-, and three-stage designs are covered as well as allocations for stratified sampling and probability proportional to size sampling. Sample allocations can be computed that minimize the variance of an estimator subject to a budget constraint or that minimize cost subject to a precision constraint. The package also contains some specialized functions for estimating variance components and design effects. Several finite populations are included that are useful for classroom instruction.


Introduction
Samples from finite populations are one of the mainstays of research in demographics, economics, and public health.In the U.S., for example, the Consumer Price Index is based on samples of business establishments and households (Bureau of Labor Statistics, 2013, chap. 17); the unemployment rate is estimated from the Current Population Survey, which is a sample of households (Bureau of Labor

Simple random sampling
Simple random sampling without replacement (srswor) is a method of probability sampling in which all samples of a given size n have the same probability of selection.The function sample in the base package in R can be used to select simple random samples either with or without replacement.One way of determining an srswor sample size is to specify that a population value θ be estimated with a certain coefficient of variation (CV) which is defined as the ratio of the standard error of the estimator, θ, to the value of the parameter: CV( θ) = Var( θ)/θ.For example, suppose that y k is a value associated with element k, U denotes the set of all elements in the universe, N is the number of elements in the population, and the population parameter to be estimated is the mean, ȳU = ∑ k∈U y k /N.With a simple random sample, this can be estimated by the sample mean, ȳs = ∑ k∈s y k /n, where s is the set The R Journal Vol.7/2, December 2015 ISSN 2073-4859 of sample elements and n is the sample size.Setting the required CV of ȳs to some desired value CV 0 in an srswor leads to a sample size of where S 2 U is the population variance of the y k 's.The term S 2 U / ȳ2 U is referred to as the unit relvariance.If y k is a 0/1 variable identifying whether an element has a characteristic or not, then S 2 U = N(N − 1) −1 p U (1 − p U ) where p U is the proportion in the population with the characteristic.The function nCont in PracTools will make the computation in (1).In a real application, the population values in (1) must be estimated from a sample or guessed based on prior knowledge.
Another way of determining a sample size is to set a tolerance for how close the estimate should be to the population value.If the tolerance (sometimes called the margin of error) is e 0 and the goal is to be within e 0 of the population mean with probability 1 − α, this translates to requiring Pr (| ȳs − ȳU | ≤ e 0 ) = 1 − α.This is equivalent to setting the half-width of a 100(1 − α)% normal approximation, two-sided confidence interval (CI) to e 0 = z 1−α/2 V ( ȳs ).The notation z denotes the 100 percentile of the standard normal distribution.The sample size required to accomplish this is One could also require that the relative absolute error, |( ȳs − ȳU )/ ȳU |, be less than e 0 with a speci- fied probability.In that case, (2) is modified by replacing S 2 U with the unit relvariance, S 2 U / ȳ2 U .Both calculations can be made using the function nContMoe in PracTools.When estimating a proportion, there are options other than a normal approximation confidence interval on p U for setting a margin of error.Two are to work with the log-odds, p U /(1 − p U ), or to use the method due to Wilson (1927), which are both available in PracTools.
Another estimand in a survey might be the difference in means or proportions.The difference could be between two disjoint groups or between the estimates for the same group at two different time periods.The standard approach in such a case would be to find a sample size that will yield a specified power for detecting a particular size of the difference.The functions power.t.test and power.prop.test in the R package stats will do this for independent simple random samples.
The case of partially overlapping samples can also be handled (e.g., see Woodward 1992).For example, persons may be surveyed at some baseline date and then followed-up at a later time.An estimate of the difference in population means may be desired, but the samples do not overlap completely because of dropouts, planned sample rotation, or nonresponse.Such non-overlap would be common in panel surveys.Suppose that s 1 and s 2 are the sets of sample units with data collected only at times 1 and 2, and that s 12 denotes the overlap.Thus, the full samples at times 1 and 2 are s 1 ∪ s 12 and s 2 ∪ s 12 .Also, suppose that the samples at the two time periods are simple random samples.Assume that the samples at times 1 and 2 are not necessarily the same size, so that n 1 = rn 2 for some positive number r.The samples might be of different sizes because of other survey goals or because the budget for data collection is different for the two times.A case that is covered by the analysis below is one where an initial sample of size n 1 is selected, a portion of these respond at time 2, and additional units are selected to obtain a total sample of size n 2 for time 2. Taking the case of simple random sampling, the difference in means at the first and second time points can be written as The variance can be expressed as where σ 2 x and σ 2 y are the population variances at the two time periods, σ xy is the element-level covariance, and n 12 is the number of units in s 12 .Writing n 12 = γn 1 and r = n 1 /n 2 , the variance becomes Var d = 1 n 1 σ 2 x + rσ 2 y − 2γrσ xy .For a one-sided test of H 0 : µ D = 0 versus H A : µ D = δ to be done with power β, the required sample size n 1 is The value of n 2 is then determined from r = n 1 /n 2 .The function nDep2sam in PracTools will perform The R Journal Vol.7/2, December 2015 ISSN 2073-4859 this calculation.nProp2sam will do a similar calculation for testing the difference in proportions with overlapping samples.

Probability proportional to size sampling
Probability proportional to size (pps) sampling can be very efficient if the measure of size (mos) used for sampling is correlated with the quantities measured in a survey.For example, enrollment in an elementary school may be related to the number of children receiving a government assistance program.In an establishment survey, the total number of employees is often correlated with other employment counts in the establishment, like the number who participate in a retirement plan.In a hospital survey, the number of inpatient beds is usually related to numbers of patients discharged in a month's time.Household samples are often selected using several stages, the first of which is a sample of geographic areas.An effective mos is typically the number of persons or housing units in each geographic area.
If a variable follows a certain regression structure in the population, then an optimal measure of size can be estimated.The key finding is due to Godambe and Joshi (1965).Isaki and Fuller (1982) extended this to a linear model M where E M (y i ) = x T i β and Var M (y i ) = v i with x i defined as a vector of x's (auxiliary variables), and β is a vector of regression slopes of the same dimension as x i .Assume that a population total is estimated and a regression estimator is used that is approximately unbiased when averaging over the model and a probability sampling design.In that case, √ v i is the best mos for pps sampling.
A model that may fit some establishment or institutional populations reasonably well has a variance with the form, Var M (y i ) = σ 2 x γ i , where x i is a mos and γ is a power.Typical values of γ are in the interval [0, 2].The function gammaFit in PracTools returns an estimate of γ using an iterative algorithm.The algorithm is based on initially running an ordinary least squares (OLS) regression of y i on x i .The OLS residuals, e i , are then used to regress log e 2 i on log (x i ) with an intercept.The slope in this regression is an estimate of γ.This procedure iterates by using the latest estimate of γ in a weighted least squares regression of y i on x i .The parameter γ is then re-estimated in the logarithmic regression.The algorithm proceeds until some user-controllable convergence criteria are met.
The variance formulas for pps without replacement sampling are difficult or impossible to use for sample size determination because they involve joint selection probabilities of units and the sample size is not readily accessible.One practical approach is to use a variance formula appropriate for pps with replacement (ppswr) sampling.The simplest estimator of the mean that is usually studied with ppswr sampling is called "p-expanded with replacement" (pwr) (Särndal et al., 1992, chap.2) and is defined as ŷpwr = 1 where p i is the probability that element i would be selected in a sample of size 1.A unit is included in the sum as many times as it is sampled.The variance of ŷpwr in ppswr sampling is where t U is the population total of y.If the desired coefficient of variation is CV 0 , Equation ( 6) can be solved to give the sample size as We later give an example of how (7) may be evaluated using PracTools.

Stratified sampling
Stratified sampling is a useful way of restricting the dispersion of a sample across groups in a population.It can also lead to improvements in precision of overall estimates if an efficient allocation to the strata is used.For example, establishments can be stratified by type of business (retail, wholesale, manufacturing, etc.).Other methods of creating strata are provided by the R package stratification (Baillargeon and Rivest, 2014).Given that strata have been created, there are various ways of efficiently allocating a sample to strata: (i) minimize the variance of an estimator given a fixed total sample size (Neyman allocation), (ii) minimize the variance of an estimator for a fixed total budget, (iii) minimize the total cost for a target CV or variance of an estimator, or (iv) allocate the sample subject to several CV or cost criteria subject to a set of constraints on stratum sample sizes or other desiderata.The last is referred to as multicriteria optimization.
The R Journal Vol.7/2, December 2015 ISSN 2073-4859 The standard texts noted earlier give closed form solutions for (i) and (ii).For example, suppose a mean is estimated, an srswor is to be selected within each stratum (h = 1, ..., H), and that the total cost can be written as C = C 0 + ∑ H h=1 c h n h where C 0 denotes fixed costs that do not vary with the sample size, c h is the cost per-element in stratum h, and n h is the number of elements sampled from stratum h.The allocation to strata that minimizes the variance of the estimated mean subject to a fixed total budget C is where W h is the proportion in the population in stratum h and S h is the population standard deviation in stratum h of the variable whose mean is estimated.This and the allocations for (i) and (iii) above can be found using strAlloc in PracTools.
Multicriteria optimization and allocations with constraints are more realistic for multipurpose surveys.In some cases, solutions to particular allocation problems are available as in Gabler et al. (2012).More generally, the alabama (Varadhan, 2015) and Rsolnp (Ghalanos and Theussl, 2014) packages will solve nonlinear optimization problems with constraints and can be very useful for complicated sample allocations.Among the constraints that are used in practical work are ones on minimum and maximum stratum sample sizes and relvariances of overall and individual stratum estimates.Theussl and Borchers (2015) present a CRAN Task View on optimization and mathematical programming in R.An R package that will form strata for multipurpose samples is SamplingStrata (Barcaroli, 2014).Also, the Solver add-on to Microsoft Excel (Fylstra et al., 1998) will handle allocation problems that are quite complex and is easy to use.Use of these tools for sample allocation is covered in some detail in Valliant et al. (2013, chap. 5).

Two-and three-stage sampling
Two-and three-stage sampling is commonplace in household surveys but can be also used in other situations.For example, the U.S. National Compensation Survey selects a three-stage samplegeographic areas, establishments, and occupations-to collect compensation data (Bureau of Labor Statistics, 2013, chap.8).Allocating the sample efficiently requires estimates of the contribution to the variance of an estimate by each stage of sampling.
As an example, consider a two-stage design in which the primary sampling units (PSUs) are selected with varying probabilities and with replacement and elements are selected at the second-stage by srswor.As noted earlier, determining sample sizes as if the PSUs are selected with-replacement is a standard workaround in applied sampling to deal with the fact that without-replacement variance formulas for pps samples are too complex to use for finding allocations.(Other selection methods along with three-stage designs are reviewed in Valliant et al. (2013, chap.9).)Let m be the number of sample PSUs, N i be the number of elements in the population for PSU i, and suppose that the same number of elements, n, is selected from each PSU.The pwr-estimator of a total is where ti = N i n ∑ k∈s i y ik is the estimated total for PSU i from a simple random sample and p i is the 1-draw selection probability of PSU i, i.e., the probability in a sample of size one.The variance of tpwr is , t i is the total of y for PSU i, and U2i is the population variance of y within PSU i. Dividing this by t 2 U and assuming that the within-PSU sampling fraction, n/N i , is negligible, we obtain the relative variance (relvariance) of tpwr as, approximately, V tpwr U2i /p i , and δ = B 2 /(B 2 + W 2 ).A simple cost function for two-stage sampling assumes that there is a cost per sample PSU and a cost per sample element of collecting and processing data.We model the total cost as The R Journal Vol.7/2, December 2015 ISSN 2073-4859 where C 0 = costs that do not depend on the number of sample PSUs and elements; C 1 = cost per sample PSU; and C 2 = cost per element within PSU.
The optimal number of units to select per PSU, i.e., the number that minimizes the approximate relvariance, is Only the ratio of the unit costs needs to be known in order to compute nopt .To find the optimal m for a fixed total cost, we substitute nopt into the cost function to obtain Alternatively, to find the optimal m for a fixed relvariance, CV 2 0 , nopt is substituted into the relvariance formula (9).clusOpt2 in PracTools will do either of these calculations.
For three-stage sampling, suppose that m PSUs are selected with varying probabilities and withreplacement, n secondary sampling units (SSUs) are selected within each PSU by srswor, and q elements are sampled by srswor within each sample SSU.This design is referred to as ppswr/srs/srs below.The relvariance of the pwr-estimator of a total in such a three-stage sample (with a negligible sampling fraction in the second and third stages) can be written as, e.g., see Hansen et al. (1953) and Valliant et al. (2013, chap.9): where The variance components and other terms in Equation ( 12) are defined as: U is the unit relvariance of y in the population with Q being the total number of elements; , the variance among the M PSU totals; is the unit variance of the N i SSU totals in PSU i with t ij = ∑ k∈U ij y k being the population total for PSU/SSU ij, tUi = ∑ j∈U i t ij N i is the average total per SSU in PSU i; U3i is the elementlevel variance among all elements in PSU i; and The terms δ 1 and δ 2 are referred to as measures of homogeneity, as is δ for two-stage sampling.Equation ( 12) is useful for sample allocation because the measures of homogeneity are in [0, 1], k 1 and k 2 are usually near 1, and Ṽ can usually be estimated.
To arrive at an optimal allocation, costs need to be considered.A cost function for three-stage sampling, analogous to the one for two-stage sampling, is The R Journal Vol.7/2, December 2015 ISSN 2073-4859 The term C 0 is again costs that do not depend on the sample sizes at different stages; C 1 is the cost per PSU; C 2 is the cost per SSU; and C 3 is the cost per element within each SSU.Minimizing the ppswr/srs/srs relvariance in Equation ( 12) subject to a fixed total cost gives the following optima: If a target relvariance is set at CV 2 0 , then the equations for finding the optima qopt and nopt are the same.The optimum number of PSUs is found by substituting nopt and qopt into the relvariance in Equation ( 12).clusOpt3 will do these computations for three-stage samples.

Two-phase sampling
In finite population sampling, a distinction is drawn between multistage sampling and multiphase sampling.In a multiphase sample, an initial sample is selected, some characteristics of the units are observed, and a decision is made about how to select a subsample from the initial sample based on what has been observed.In multistage sampling, the same design is used in later stages regardless of what was found in the first-stage units.There is a more technical definition of the difference between multistage and multiphase, but it is unimportant for this discussion.An example of two-phase sampling is to select a subsample of nonrespondents to the initial phase to attempt to get them to cooperate.This is known as a nonresponse follow-up study (NRFU).
Another type of two-phase design is double sampling for stratification.In this design, information is collected in the first phase which is then used to stratify elements for second phase sampling.For example, researchers working to develop a case definition for undiagnosed medical symptoms in U.S. personnel serving in the 1991 Persian Gulf War surveyed a stratified simple random sample of Gulf War-era veterans (Iannacchione et al., 2011).Based on survey responses in the first phase, respondents were classified as likely having or not having a certain type of illness.Blood specimens were requested from randomly sampled phase-1 respondents within the illness strata and analyzed using expensive tests.
As an illustration, take the case of double sampling for stratification.Cochran (1977) and Neyman (1938) give the two-phase variance of an estimated mean or proportion when phase-1 is a simple random sample of n (1) elements, phase-2 is a stratified simple random sample (stsrs) of n (2) elements, and an optimal allocation to strata is used in the second phase.The sampling fractions at both stages are assumed to be negligible.The optimal proportion of the phase-2 sample to assign to stratum h for estimating the population mean is n (2)h /n (2) = W h S h / ∑ h W h S h .The formula for the variance of an estimated stratified mean with this allocation is where ȳUh is the stratum h population mean.The phase-2 subsampling rate from the phase-1 sample units that minimizes Equation ( 17) is where c (1) and c (2) are the per-unit costs in the first and second-phases, respectively.The formulas for the phase-1 and phase-2 sample sizes that minimize V opt subject to a fixed total cost C are The function dub in PracTools will calculate the optimal second-phase sampling fraction and the phase-1 and phase-2 sample sizes that will minimize the variance in Equation ( 17).The function also computes the size of an srs that would cost the same as the two-phase sample and the ratio of the The R Journal Vol.7/2, December 2015 ISSN 2073-4859 two-phase stratified variance to the srs variance.

Extension to nonlinear estimators and limitations
PracTools will calculate sample sizes for estimators whose variance can be written in one of the forms given in Section Designing survey samples.This directly covers linear estimators of means and totals.
A nonlinear estimator that is a differentiable function of a vector of estimated totals is also covered.But, a user must do some work to linearize the estimator and determine the inputs that are required for a PracTools function.Suppose that the estimator is θ = f t1 , ..., tp where f is a differentiable function and tj is a linear estimator of a population total, t j (1, 2, ..., p).The linear approximation to θ − θ is where θ = f t 1 , ..., t p , and the partials are evaluated at the population values.The repeated sampling variance of this approximation is then the same as the variance of ∑ p j=1 ∂ f ∂t j tj since θ and t j are treated as constants.Taking the case of two-stage sampling, suppose that the estimator of the total for variable j is tj = ∑ i∈s ∑ k∈s i w ik y ik (j) where w ik is a weight for element k in PSU i, y ik (j) is its data value, s is the set of sample PSUs, and s i is the sample of elements within SSU i. Substituting this into (18) and reversing the order of summation between PSUs and variables leads to the expression where ti (j) = ∑ k∈s i w ik y ik (j).One then computes the design-variance of ẑ based on the particular sample design used and the form of the derivatives.This general approach is known as the linear substitute method and is described in detail in Wolter (2007).
In a two-stage sample where m PSUs are selected with replacement and with varying probabilities, and n elements are selected by simple random sampling from each PSU, (9) applies.The weight is defined as w ik = (mp i ) −1 (N i / n) where p i is the 1-draw selection probability, as before.If the estimator is the mean computed as θ = ty / M with M = ∑ i∈s ∑ k∈s i w ik , then ẑi = (Mmp i ) −1 N i e i where e i = ȳi − ȳU with y i being the sample mean of y in PSU i and ȳU the population mean.
with tzi = N i / n ∑ k∈s i e ik and e ik = y ik − ȳU .Expression (9) would then be evaluated with e ik replacing y ik .With the linear substitute method, residuals typically appear in the linear approximation and are the basis for a variance estimator.Population quantities, like M and ȳU , are replaced by sample estimates in a variance estimator.The Section Examples gives an illustration of this method using one of the datasets in PracTools.
Nonetheless, there are types of estimators that our package does not cover.Quantile estimators require a special approximation and variance formula (Francisco and Fuller, 1991) that does not come from the standard linearization approach.The Gini coefficient, used as a measure of income inequality (e.g., Deaton, 1997), is another example of an estimator that is too complicated to linearize using the methods above.

The R package PracTools
PracTools is a collection of specialized functions written in R along with several example finite populations that can be used for teaching.The package is available from the Comprehensive R Archive Network (CRAN) at http://CRAN.R-project.org/package=PracTools.Because the function code is visible, the routines can be modified (and improved) by any user.A brief description of the functions is given in Table 1.The use of the functions, their input parameters, and the values they return are described in the help files.The package also contains nine example populations of different types.These are listed in Table 2.

BW2stagePPS, BW2stageSRS
Variance components in two-stage samples from a population frame BW2stagePPSe Estimated variance components in two-stage samples from a sample BW3stagePPS Variance components in three-stage samples from a population frame BW3stagePPSe Estimated variance components in three-stage samples from a sample clusOpt2 Optimal allocation in a two-stage sample clusOpt2fixedPSU Optimal second-stage sample size in a two-stage sample when the PSU sample is fixed clusOpt3 Optimal allocation in a three-stage sample clusOpt3fixedPSU Optimal second and third-stage sample sizes in a three-stage sample when the PSU sample is fixed CVcalc2, CVcalc3 Compute the coefficient of variation of an estimated total in two-and three-stage designs deffH, deffK, deffS

Examples
This section gives some examples for computing sample sizes for estimating proportions and differences of proportions, an allocation to strata, and the optimal numbers of PSUs, secondary units, and elements in a three-stage sample.
The R Journal Vol.7/2, December 2015 ISSN 2073-4859 nDep2sam will also accept a vector of differences as input, e.g., del <-seq(2,10,0.001).This makes it easy to generate a plot like in Figure 1, where the sample size in each group is plotted against the difference in means for several levels of power.This is a useful way to present options to users.

Probability proportional to size sampling
The function gammaFit will estimate the variance parameter in the model The code below uses the hospital population bundled with PracTools to estimate the variance power in the model, E M (y) = β 1 √ x + β 2 x, V M (y) = σ 2 x γ , returning γ = 1.612009.Using the function UPrandomsystematic from the package sampling, a sample of size 30 is then selected with probability proportional to √ x γ, which is the optimal measure of size for estimating the population total of y.

Allocations in two-and three-stage samples
When designing multistage samples, decisions must be made about how many units to select at each stage.To illustrate this, we consider a three-stage sample.A considerable amount of data is needed to estimate realistic ingredients required for clusOpt3.The function, BW3stagePPSe, will estimate B 2 , W 2 , W 2 2 , W 2 3 , δ 1 , and δ 2 from a three-stage where the first-stage is selected ppswr and the last two stages are selected by srswor.BW2stagePPSe does similar calculations for two-stage sampling.Variance component estimation is, of course, a difficult area where a number of alternatives have been developed in the model-based literature.The forms used in BW2stagePPSe and BW3stagePPSe are fairly simple ANOVA-type estimates.These estimates have known defects, like occasionally being negative.
The following example computes a three-stage allocation that minimizes the variance of the pwrestimator assuming that the budget for variable cost is 100,000; the PSU, SSU, and per-element costs are 500, 100, and 120, respectively; δ 1 = 0.01, δ 2 = 0.10; the unit relvariance is Ṽ = 1; and the ratios, k 1 and k 2 , are both 1. cal.sw = 1 specifies that the optima be found for a fixed total budget.The full description of the input parameters can be found in the help file for clusOpt3.> clusOpt3(unit.cost= c(500, 100, 120), delta1 = 0.01, delta2 = 0.10, unit.rv= 1, + k1 = 1, k2 = 1, tot.cost = 100000, cal.sw = 1) C1 = 500 C2 = 100 C3 = 120 delta1 = 0.01 delta2 = 0.1 unit relvar = 1 k1 = 1 k2 = 1 cost = 1e+05 m.opt = 28.3n.opt = 7.1 q.opt = 2.7 CV = 0.0499 Along with the inputs, the output includes the optimal sample size for each stage and the CV that is anticipated for the pwr-estimator given that design.The sample sizes above can be rounded and then used in the sampling package to select units at each of the stages.For example, suppose that 28 PSUs and 7 SSUs will be selected with probabilities proportional to a mos.Within each sample SSU, a sample of 3 elements will be selected via srswor.The PSUs can be selected using the function cluster and the data extracted.Then, a cluster sample of 7 SSUs can be selected from each of those 28 units in a loop, again using cluster and the data for those sample SSUs extracted.The sample of SSUs would then be treated as strata and the strata function used to select 3 elements from each SSU using srswor.

Double sampling for stratification
The function dub will compute the allocation to strata for a double sampling design in which phase-1 is used to assign units to strata.The function takes these input parameters: c1 cost per unit in phase-1 c2 cost per unit in phase-2 Ctot total variable cost Nh vector of stratum population counts or proportions Sh vector of stratum population standard deviations Yh.bar vector of stratum population means The inputs, N h , S h , and Ȳh , will typically have to be estimated from the first-phase sample.The example below computes the allocation to four strata assuming a total cost of 20,000 and unit costs of c (1) = 10 and c (2) = 50.A proportion is being estimated.

Sample size for a nonlinear estimator
To illustrate a calculation for a nonlinear estimator, consider the proportion of Hispanics with insurance coverage in the MDarea.pop,which is part of the package.Define y 2k to be 1 if a person is Hispanic and 0 if not; α 1k = 1 if a person has insurance coverage.Then, y 1k = α 1k y 2k is 1 if person k has insurance and is Hispanic and is zero otherwise.The linear substitute is z k = y 1k − θy 2k where θ is the proportion of Hispanics with insurance coverage.In this case, z k can take only three values: −θ, 0, and 1 − θ.If a simple random sample of clusters and persons within clusters is selected, BW2stageSRS can be used to compute B 2 , W 2 , and δ using the linear substitutes as inputs.Assuming that the full population is available, the R code is the following.We do the calculation for clusters defined as tracts (a small geographic area with about 4000 persons defined for census-taking).

Figure 1 :
Figure 1: Sample size in each group for detecting a range of differences in means with several levels of power.