{'fraction_non_alphanumeric': 0.031243159340112635, 'fraction_numerical': 0.009969950846749318, 'mean_word_length': 4.607856266041736, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Inference under Cluster Sampling with Probability Proportional to Size

Cluster sampling is common in survey practice, and the corresponding inference has been predominantly design-based. We develop a Bayesian framework for cluster sampling and account for the design effect in the outcome modeling. We consider a two-stage cluster sampling design where the clusters are first selected with probability proportional to cluster size, and then units are randomly sampled inside selected clusters. Challenges arise when the sizes of nonsampled cluster are unknown. We propose nonparametric and parametric Bayesian approaches for predicting the unknown cluster sizes, with this inference performed simultaneously with the model for survey outcome. Simulation studies show that the integrated Bayesian approach outperforms classical methods with efficiency gains. We use Stan for computing and apply the proposal to the Fragile Families and Child Wellbeing study as an illustration of complex survey inference in health surveys.


Introduction
Cluster sampling has been widely implemented in epidemiology and public health surveys (Carlin and Hocking, 1999). We develop a Bayesian approach for survey inference under cluster sampling, particularly in the absence of design information for nonsampled clusters. The potential advantages of Bayesian methods are for small area estimation and adjusting for many poststratification factors (Gelman, 2007). However, most of the work in this area has been done for one-stage sampling or ignoring clustering in the data collection; thus it would be a potentially useful contribution to bring cluster sampling into the class of problems for which there are available Bayesian models. In the present paper, we consider two-stage probability proportional to size (PPS) sampling, with the understanding that other designs could be modeled in similar ways.
Cluster sampling increases cost efficiency when partial clusters are included in the probability sampling framework. Bayesian cluster sampling inference is essentially outcome prediction for nonsampled units in the sampled clusters and all units in the nonsampled clusters. The design information should be accounted for in the modeling, but design information for nonsampled clusters is often unknown or inaccessible. We introduce estimation strategies for such information and connect multilevel regression models to sampling design as a unified framework for survey inference.
We consider the design that involves first sampling primary sampling units (PSUs) and then sampling secondary sampling units (SSUs) within selected PSUs. The two-stage cluster sampling design has played an important role in the sequential data collection process of many big health surveys, such as the National Health Interview Survey and the Medical Expenditure Panel Survey. This design requires a complete listing of PSUs and a complete listing of SSUs only within selected PSUs, and thus is widely used when generating a sampling frame of every unit in the population is infeasible or impractical. For example, in designing a nationally representative household survey, generating a complete listing of every household in the country requires essentially as much effort as a complete census of all households. Instead, the sampling proceeds in stages, first sampling PSUs such as counties, cities, or census tracts. The PSUs are sampled with probability proportional to size, which is commonly the number of SSUs in the PSU but can be a more general measure of size, such as annual revenue or agricultural yield. SSUs are then randomly selected within selected PSUs, often with a fixed number or proportion. This design assumes independence and invariance of the second-stage sampling design . Invariance means that the sampling of SSUs is independent of which PSUs are sampled, and independence means sampling of SSUs in one PSU is independent of sampling in other PSUs. For clarification, a two-phase design is one in which one or both assumptions do not hold.
Our motivating application survey, the Fragile Families and Child Wellbeing study (Reichmann et al., 2001), was collected via a multistage design, where two-stage cluster sampling served as a key step. The study aims to examine the conditions and capabilities of new unwed parents, the wellbeing of their children and the policy and environmental effect. To obtain a nationally representative sample of non-marital births in large U.S. cities, the study sequentially sampled cities, hospitals, and births. The sampling of cities used a stratified random sample of all U.S. cities with 200,000 or more people, where the stratification was based on policy environments and labor market conditions in the different cities. Inside each stratum, cities were selected with probability proportional to the city population size. In the selected cities, all hospitals in the small cities were included, while a random sample of hospitals or the hospital with the largest number of non-marital births was selected in large cities. Lastly, a predetermined number of births were selected inside each hospital. Classical weighting adjustment for the complex study design results in highly variable weights (Carlson, 2008), thus yielding unstable inferences. We would like to develop a Bayesian framework to account for the complex designs under two-stage cluster sampling.
Our goal is to develop hierarchical models and account for design effects to yield valid and robust survey inference. Bayesian hierarchical models are well-equipped to handle the multistage design and stabilize estimation via smoothing. As an intermediate step, twostage cluster sampling is crucial in the Fragile Families study to select cities and hospitals. However, cluster sampling presents unique methodology challenges as little information is available on the unselected clusters. This article uses the Fragile Families study as an illustration and focuses on Bayesian cluster sampling inference to build a unified survey inference framework. The unified framework can be extended under a complex sampling design, as discussed in Section 5.
We illustrate finite population inference with the estimation of the population mean in a two-stage cluster sample. Specifically, we consider a population of J clusters, with each cluster j containing N j units and a total population size of N = J j=1 N j . Let I j denote the inclusion indicator for cluster j and I i|j denote the inclusion indicator for unit i in cluster j, i = 1, . . . , N j [i] , where j[i] denotes the cluster to which unit i belongs. Clusters are sampled with probability proportional to the measure of size M j , which is known to the analyst only for the sampled clusters. Our goal is to estimate the finite population mean of the survey variable y, which, for a continuous variable is defined as where y j represents the mean of y in cluster j. For a binary outcome, we seek to estimate the population proportion, which is given by where y (j) is the population total in cluster j. Classically, inference in survey sampling has been design-based. The design-based approach treats the survey outcome y as fixed, with randomness arising solely from the random distribution of the inclusion indicator I. Design-based estimators have the advantage of being design-consistent, where design-consistency means that the estimator will converge to the true value as the population and sample sizes increase under the given sampling design. However, they are often unstable with large standard errors. For estimating the finite population mean of an outcome y i , the classical design-based estimator for a single-stage sample s of size n is the Hájek estimator where π i is the inclusion (selection and response) probability of unit i. In the two-stage sample s, when J s out of J clusters are selected with n j 's sampled SSUs, for convenience labeled as j = 1, . . . , J s , the estimator becomes where π j is the selection probability of cluster j, and π i|j is the selection probability of unit i in cluster j given that cluster j was sampled. The design-based approach does not require a statistical model for the survey outcomes but implicitly assumes the specific outcome model structure and linearity (Chen et al., 2017). The performance relies on the validity of the model assumptions, and then the estimators can yield biased inference under invalid assumptions. Another major challenge with design-based estimators comes in estimating their variance. The variance of a design-based estimator generally requires knowledge of not only the inclusion probability π i for a given unit i, but also the joint inclusion probability π ii for any two units i and i in the population. This information is often unavailable in practice, such as the unknown measure of size for nonsampled clusters under the PPS setting. Joint inclusion probabilities can be challenging to compute even for straightforward sampling designs, and variance estimators for designbased estimators are often based on simplifications and approximations. Furthermore, the inverse-probability of inclusion weighting often leads to highly variable estimators.
Bayesian inference, in contrast, directly models both the inclusion indicators I i and the survey outcomes y i . The Bayesian approach to survey inference has many advantages over the design-based approach, including the ability to handle complex design features such as multistage clustering and stratification, stabilized inference for small-sample problems, incorporation of prior information, and large-sample efficiency (Little, 2004). When the design variables are included in the model, the selection mechanism becomes ignorable (Rubin, 1983;Gelman et al., 2013), and we can model the outcomes y alone, instead of jointly modeling y and the inclusion indicator I. The importance of including design variables in the model has also been emphasized for missing data imputation (Schafer, 1997;Reiter et al., 2006).
Unfortunately, in many (arguably most) practical situations, the set of design variables is not available for the entire population and is instead known only for sampled clusters or units. In the case of PPS sampling, where the design variables consist of the cluster measures of size {M j } J j=1 , we as the survey analysts may only have access to M j (or, equivalently, the inclusion probability π j ) for the sampled subset of J s clusters. This missing information on measures of size causes methodology challenge in the Bayesian setting because we cannot predict the values of y for the nonsampled clusters without it. We need to model the values of M j for nonsampled clusters before we are able to make inferences about y conditional on the design information.
Existing Bayesian approaches to this problem (Zangeneh et al., 2011;Zangeneh and Little, 2015) consider the case of single-stage PPS sampling. In addition, they separate estimation of the missing measure sizes and inference for the finite population quantities into two steps. In contrast, we propose an approach that integrates these steps into one model for a two-stage cluster sample. Our model allows for both cluster-and unit-level information to be used when both are available in certain cases. For much of this paper, we assume the measure of size is equal to the cluster size N j and use N j in place of M j for simple illustration.
The rest of this paper proceeds as follows. Section 2 first gives an overview of current approaches to finite population inference under PPS and then describes our approach and its advantages. In Section 3, we describe a simulation study to investigate the performance of our method and compare with other literature methods. We apply our proposal to the Fragile Families study in Section 4 and discuss the conclusions and extensions in Section 5.

Methods
In the two-stage cluster sampling, a fixed number J s of clusters are sampled with PPS, so that the probability of cluster j being included in the sample is proportional to N j : We only observe N j 's for the clusters in the sample, that is, the empirical distribution of (N j |I j = 1). Our proposed procedure simultaneously models the population cluster sizes and the outcome and propagates the estimation uncertainty.
Let x i denote the auxiliary variables that are predictive for the outcome. The observed data are (y obs , x obs , N obs , x 1:J , N, J, J s ), where x 1:J is the cluster-level mean of the covariate x for all clusters j = 1, . . . , J, and N , J, and J s are the total population size, total number of clusters, and number of sampled clusters, respectively. The subscript obs denotes the observed portions of the variables: y obs = {y i : i = 1, . . . , n j[i] , j = 1, . . . , J s }, x obs = {x i : i = 1, . . . , n j[i] , j = 1, . . . , J s }, N obs = {N j : j = 1, . . . , J s }, where for convenience we number the sampled clusters j = 1, . . . , J s and the nonsampled clusters as j = J s + 1, . . . , J. We assume that x i is known for all sampled units, and that x j is known for all clusters. If x is a demographic covariate, in practice it's often the case that we know demographic characteristics of clusters even if the cluster size is unknown.
The goal is to estimate the finite population mean y, defined for a continuous outcome, where y obs,j is the mean of the sampled units in sampled cluster j, y exc,j is the mean of the nonsampled units in cluster j, and N exc,j is the size of nonsampled cluster j that is unknown. For a binary outcome, the population proportion is where y (j) is the total of all units in cluster j, y obs,(j) is the total of sampled units in sampled cluster j and y exc,(j) is the total of the binary outcome in nonsampled units in cluster j. We assume the continuous survey outcome y is related to the covariate x and cluster sizes N j in the following way: where φ are the parameters governing the distribution of the cluster sizes N j . The model assumes the regression coefficients are cluster-varying and depend on the cluster sizes. We use random-effects model to borrow information across clusters. While fixed-effects model with cluster membership indicators can also be used to quantify the cluster effect, fixed cluster effects models may increase the variance, as shown by Reiter et al. (2006) and Andridge (2011) in the context of missing data imputation. In addition, predictions cannot be made for nonsampled clusters using fixed-effects models since no units are available.
Our model for a binary outcome is identical, except that we modify (4) to be and omit (6). We exclude unit-level covariates in the binary case because the nonlinear nature of the inverse logit link function makes it challenging to make use of data at the unit level. Specifically, predicting y exc,j requires knowing x i for all nonsampled units in cluster j, and if we knew this, clearly we would have also known N j for nonsampled clusters j. We use the centered logarithms of the cluster sizes log c (N j ) as predictors; we work on the logarithm scale to better accommodate large cluster sizes and center for interpretation convenience. The sampling is assumed to be ignorable after including the design variables in the outcome model.
We assign an estimation model p(N j | φ) to the cluster sizes, which we observe only for the sampled clusters. We develop both nonparametric and parametric modeling strategies to predict the cluster sizes of nonsampled clusters.
Because of the independence, invariance and ignorability assumptions in the two-stage cluster sampling, the distribution of the outcome y, given the design variables, is the same in the sample and the population; that is, the observed data likelihood is the same as the complete data likelihood, p(y obs | x obs , N obs , ψ) = p(y | x, N, I = 1, ψ) = p(y | x, N, ψ).
The challenge lies in estimating the distribution of the N j 's when the sampling is informative. Under PPS sampling, the probability of observing a cluster of size N j is We consider both nonparametric and parametric modeling strategies for the prior distribution p(N j ) (also called the population distribution, to distinguish from the distributions of sampled and nonsampled cluster sizes) in (7). First, we introduce the Bayesian bootstrap algorithm in Section 2.1 as a nonparametric approach to predicting the unobserved N j 's. Second, we investigate two parametric distributional assumptions in Section 2.2 for p(N j ), the negative binomial and lognormal distributions. Here our goal is to directly model the distribution of the cluster sizes accounting for the fact that the observed distribution is biased from the complete population distribution. Following Patil and Rao (1978), we refer to these parametric choices as size-biased distributions.

Bayesian bootstrap
For a nonparametric model of the sampled cluster sizes, we implement the Bayesian bootstrap algorithm in Little and Zheng (2007) and Zangeneh et al. (2011) for one-stage PPS sampling and modify it under two-stage PPS sampling. Without a parametric assumption for p(N j ), we connect p(N j | I j = 0) with p(N j | I j = 1) through the empirical distributions under PPS sampling. Assume the N j 's observed for the sampled clusters have B unique values N * 1 , . . . , N * B , and let k 1 , . . . , k B be the corresponding counts of these unique sizes, such that b k b = J s . Let ψ b denote the probability of observing a cluster of size N * b in the sample: . We can then model the counts k = (k 1 , . . . , k B ) as multinomially distributed with total J s and parameters ψ = (ψ 1 , . . . , ψ B ). The observed likelihood L obs (ψ) is, where I(·) is an indicator function, I(·) = 1 if the inside expression is true and 0 otherwise. The ψ's are given a noninformative Haldane prior: p(ψ 1 , . . . , ψ B ) = Dirichlet(0, . . . , 0), a conjugate Dirichlet prior distribution. The posterior distribution of ψ is then Suppose the unique values of N j 's cover all possible values in the population. Assume k b is the number of nonsampled clusters with size N * b , for b = 1, . . . , B, and let ψ b denote the probability of an unobserved cluster having size N . Then the counts of the B unique sizes among the nonsampled clusters, (k 1 , . . . , k B ), follow a multinomial distribution with total J − J s = b k b and probabilities (ψ 1 , . . . , ψ B ): Using Bayes' rule, we can write ψ b as is the conditional cluster selection probability known in the PPS sample, J s is the number of sampled clusters, and N is the population size. This approach essentially adjusts the probability of resampling an observed size N * b by the odds of a cluster of that size not being sampled, so that smaller sizes are upweighted relative to larger ones.
Given the posterior draws of ψ b 's and k b 's, we create k * b replicates of the size N * b , yielding a sample of the nonsampled cluster sizes from their posterior predictive distribution.
The Bayesian bootstrap for cluster sampling is similar to the "two-stage Pólya posterior" approach proposed by Meeden (1999), which simulates draws that form a population of clusters and then an entire population of elements within each cluster. Zhou et al. (2016) incorporated weights in Bayesian bootstrap for multiple imputation in two-stage cluster samples. Si et al. (2015) uses a similar approach to estimating the poststrafication cell sizes constructed by the survey weights.
The Bayesian bootstrap avoids parametric assumption on the population distribution p(N j ) and use the empirical distribution in the observed clusters. This implicitly introduces a noninformative prior distribution on N j 's. However, this approach restricts the draws for the nonsampled cluster sizes to come from the set of observed cluster sizes, where small clusters may be omitted under PPS sampling. While the Bayesian bootstrap is a robust algorithm for predicting the unknown N j 's, we can achieve efficiency gains with a parametric distribution on p(N j ), especially when prior distribution information is available.

Size-biased distributions
Inducing parametric sized-biased distributions follows the superpopulation concept in the model-based survey inference literature and incorporates informative prior information. In practice, we may have some knowledge about the cluster sizes, such as the distribution in a similar population or from previous years. We can incorporate this additional information through the prior distribution specification. Sized-biased distributions were discussed by Patil and Rao (1978) for population size estimation. We consider a discrete and a continuous distribution as candidates for modeling the size distributions. The observed likelihood is connected with the proposed population distribution via (9). Using the PPS sample, we can estimate the parameters in the population distribution and then predict the nonsampled cluster sizes.
For the discrete case, we assume the population cluster sizes N j follow a negative binomial distribution: N j ∼ NegBin(k, p), with k > 0 and p ∈ (0, 1). By normalizing the distribution in (9) and completing the algebra shown as below, we see that the sizes in the PPS sample can be written as N j = 1 + W j , where W j ∼ NegBin(k + 1, p).
For the continuous case, we use the lognormal distribution. If the population distribution is N j ∼ lognormal(µ, τ 2 ), then (N j | I j = 1) ∼ lognormal(µ + τ 2 , τ 2 ). To see this, recall that p(N j ) denotes the pdf of size variables N j in the population. Then the pdf of N j in the PPS sample is We can now simplify the denominator: Now, substitute (11) for the denominator in (10): Thus, the distribution of sampled cluster sizes in the PPS sample is (N j |I j = 1) ∼ lognormal(µ+ τ 2 , τ 2 ).
Regardless of the parametric model we choose, in order to generate predictions of the nonsampled cluster sizes, we need to draw from p(N j | I j = 0). We apply rejection sampling and use samples from p(N j ) to approximate the sampling from p(N j | I j = 0).
where G Pr(I j = 0 | N j )/Pr(I j = 0) has a constant upper bound shown as below. The marginal probability selection for cluster j is Pr(I j = 1) = J s /J, and the joint distribution of (N j , I j ) in the PPS sample is p(N j , I j = 1) = cN j p(N j ), where c is a constant. And Hence, c = J s /(JE(N j )). Then Zangeneh et al. (2011) assumes that E(N j ) = N/J, approximated by the finite sample cluster size, such that Given the posterior distribution of p(N j | −), we use rejection sampling to obtain posterior predictive samples from p(N j | I j = 0, −).

Prior specification and computation
We use the following weakly informative prior distributions as recommended by Gelman (2006), Here Cauchy + (0, 2.5) denotes a Cauchy distribution with location 0 and scale 2.5 restricted to positive values. The weakly informative prior specification will allow the group-level variance parameters to be close to 0 and have large tail values. For the parameters governing the distribution of N j , such as (k, p) in the negative binomial distribution or (µ, τ ) in the lognormal distribution, we can use noninformative priors when the number of clusters sampled is large. However, when only a few clusters are sampled, we need informative prior information to counteract the sparsity of the data and stabilize the inference. This is particularly true when using a model for the cluster sizes that includes implicit assumptions about the data. For example, as an overdispersed extension of the Poisson distribution, the negative binomial distribution assumes that the data come from a distribution whose mean is smaller than the variance. However, in a sample of only five clusters, it may well be that the sample mean is larger than the sample variance, making it difficult to fit the negative binomial distribution to the data without strong prior information. In this case, we reparameterize the negative binomial as a Gamma mixture of Poisson distributions and place a prior on the coefficient of variation (CV), the standard deviation divided by the mean. In this case, the CV works out to the reciprocal of the square root of the scale parameter of the Gamma distribution. With a small number of clusters, we expect the CV to be close to one and therefore use an exponential prior distribution with rate 1. For the lognormal distribution, we place a Cauchy + (0, 2.5) prior on the scale parameter τ . To aid estimation for the case with only a few sampled clusters, we standardize the log of the sampled cluster sizes by subtracting their mean and dividing by the standard deviation.
For the continuous outcome, in nonsampled clusters j, the posterior predictive distribution for y exc,j is (y exc,j | ·) ∼ N β 0j + β 1j x j , σ 2 y /N j , where we assume x j is known. Specifically, we draw new values of β 0j , β 1j , σ y , and N j from their posterior distributions and then draw y exc,j from the above distribution. In sampled clusters, the posterior predictive distribution for the nonsampled units is When N j is large compared to n j , as is the case in many large-scale surveys and specifically in the Fragile Families study, y exc,j is close to the cluster mean y j and is well approximated by β 0j + β 1j x j , which we calculate using the posterior means of β 0j and β 1j . The posterior computation is implemented in Stan (Stan Development Team, 2016b), which conducts full Bayesian inference and generates the posterior samples. The estimation for the outcome model and the cluster size model can be integrated into the posterior computation, which allows for uncertainty propagation throughout the parameter estimates, in contrast to previous approaches (e.g., Little and Zheng, 2007;Zangeneh and Little, 2015).
Stan is unique in providing detailed warnings and diagnostics to inform the user when posterior inferences may be unreliable due to difficulties in sampling. Divergent transitions indicate that the sampler is unable to explore a portion of the parameter space, which can lead to significant bias in the resulting posterior distribution and ultimately unreliable inferences (Stan Development Team, 2016c). Stan reports the number of divergent transitions for each run, and even one divergent transition indicates that the results may be suspect. If divergent transitions occur, we follow the recommendation of Stan developers and iteratively increase the target acceptance rate adapt delta (Stan Development Team, 2016a). If divergent transitions occur even with adapt delta = 0.99999, we switch to the noncentral parameterization and follow the same procedure for increasing adapt delta as necessary. The noncentral parameterization is a mathematically equivalent formulation for the model that can avoid posterior geometries that are difficult for HMC to explore; see Betancourt and Girolami (2013) and Stan Development Team (2016c).
To understand the importance of explicitly controlling for all design variables in this context, we also fit a model similar to (4)-(7) but with γ 0 and γ 1 set to 0. Such a model accounts for the hierarchical cluster nature of the data by allowing β 0 and β 1 to vary by cluster, but does not account for the PPS sampling design since the cluster sizes N j are excluded from the model: β 1j ∼ N(α 1 , σ 2 β 1 )

Simulation study
We perform a simulation study to compare the performance of our integrated approaches with classical design-based estimators on the statistical validity of the finite population inference. We generate a population from which we take repeated two-stage cluster samples under PPS and use each of the methods to estimate y. The population consists of J = 100 clusters, with cluster sizes N j drawn from one of two distributions. The first is a Poisson distribution with rate 500. The second is a multinomial distribution over scaled Gammadistributed sizes. Specifically, we draw J = 100 candidate cluster sizes N j as N j = 100G j , where G j ∼ Gamma(10, 1). We then take a multinomial draw from these 100 unique sizes, with the J-vector of probabilities drawn from a Dirichlet distribution with concentration parameter 10, which disperses probability mass equally across the J = 100 components. In both cases, to avoid clusters that would be selected with probability 1, we resample the J cluster sizes until none are so large to be selected with certainty. For continuous outcome, we simulate the population unit value y i from the following: where N + (µ, σ) denotes the positive part of the normal distribution with mean µ and standard deviation σ. The model for binary y is identical, except that the first line of (13) is replaced with y i ∼ Bern(logit −1 (β 0j[i] )) (and we omit β 1j ). We use the same outcome model for data generation and estimation to focus on the performance evaluation of different approaches accounting for the design effect and avoid potential model misspecification. In practice the outcome model can be adapted with flexible choices, as discussed in Section 5. We generate x i by sampling from the discrete uniform distribution between 20 and 45 and center it by subtracting the mean. We assume that x i is known for all sampled units, and that x j is known for all clusters. The cluster sizes N j 's are only known in the sampled clusters.
We sample J s < J clusters using random systematic PPS sampling with probability proportional to the cluster size N j and n j units via SRS in each selected cluster j. We consider values of J s ∈ {10, 50} and n j ∈ {0.1N j , 0.5N j , 10, 50}. Note that when n j ∈ {10, 50}, the sample is self-weighting, meaning each unit has an equal probability of selection. To see this, recall that the probability of sampling cluster j is π j ∝ N j . Since within-cluster sampling is done with SRS, the probability of sampling unit i given cluster j is selected is π i|j = n j /N j = n/N j when n j is the same for all clusters. The marginal probability of sampling unit i is therefore π i = π j π i|j ∝ N j · (n/N j ) = n, which is constant across units and clusters. Even though the final weights are constant, our studies show that the design features should be accounted in the outcome model.
For each combination of J s and n j , we draw 100 two-stage samples from the population. For each sample, we estimate the finite population mean using the methods described below.
• negbin: The negative binomial size-biased distribution as described in Section 2; • lognormal: The lognormal size-biased distribution as described in Section 2; • bb: The Bayesian bootstrap as described in Section 2; • Hájek: The Hájek estimator in (3); • greg: The generalized regression estimator (Deville and Särndal, 1992), which leverages a unit-level covariate to improve prediction. We only use this estimator for continuous y. To estimate the variances of the Hájek and generalized regression estimators, we use the formulas given in Chapter 8 of Särndal et al. (1992); 1 • cluster inds: The model in (12), which accounts for the hierarchical nature of the data via random cluster effects but does not use the cluster sizes as a cluster-level predictor in modeling β 0j and β 1j ; • knowsizes: The model in (4)- (6), where we additionally assume the cluster sizes are known for the entire population. This is the best scenario and will serve as a benchmark for the other Bayesian methods.
There are three main comparisons that we make in evaluating the results of the simulation study. First, we measure the performance of our proposed integrated Bayesian approach against that of the classical design-based estimators; we do this by comparing the performance of negbin, lognormal, and bb to that of Hájek and greg. Second, among the Bayesian methods, we want to understand when the parametric models negbin and lognormal outperform the nonparametric Bayesian bootstrap bb. Third, we compare the performances of cluster inds and knowsizes in order to understand the importance of explicitly including cluster sizes as cluster-level predictors in (5) and (6). In this case, we assume that cluster sizes are known for all clusters in the population and focus on the effects of incorrectly excluding or including the cluster sizes as cluster-level predictors in the model.
We carefully monitor the diagnostics of computation performance for each drawn sample. If divergent transitions remain, we discard the sample. We monitor the estimated potential scale reduction factor R for each parameter. This diagnostic assesses the mixing of the chains; at convergence, R = 1. If R ≥ 1.1 for any parameter, we increase the number of iterations by 1000 until all values of R are less than 1.1, up to 4000 iterations. If values of R ≥ 1.1 remain with 4000 iterations, we discard the drawn sample. The results presented here are based on a minimum of 85 simulation draws for each combination of number of clusters sampled and number of units sampled. That is, we repeatedly draw 100 samples from the population and keep the L cases with good computation performance, 85 ≤ L ≤ 100.
The results of the simulation study are in Figures 1 to 4, with each figure displaying a different combination of outcome type (continuous or binary) and population cluster size model (Poisson or multinomial). In each figure, there are six panels displaying the six metrics with which we evaluate the methods: relative bias, relative root mean squared error (RRMSE), coverage of 50% and 95% uncertainty intervals, and the average relative widths of the 50% and 95% uncertainty intervals. The relative bias is calculated as 1 L L =1 y− y y , where y is the true population mean, y is the estimated value from the -th simulation, and L is the number of simulations. RRMSE is calculated as For the Bayesian methods negbin, lognormal, bb, cluster inds and knowsizes, the 50% (95%) 1 In some cases, the sample size is so large as to make calculating the design-based variance under a non-self-weighting design difficult. This is due to the∆ k term in equations 8.6.3 and 8.9.27 in Särndal et al. (1992), which requires generating an n×n matrix, where n = Js j=1 n j . When J s = 50 and n j = 0.5N j , n can easily be 20000 or larger, making the matrix prohibitively large to compute. In these cases, we estimate the variance by randomly selecting 100 units via SRS in each sampled cluster and using those units to compute the required matrix. intervals are calculated from the 25th and 75th (2.5th and 97.5th) percentiles of the posterior predictive distribution for y. For the classical methods, we rely on asymptotic normal theory and the variance estimators given in Chapter 8 of Särndal et al. (1992). The relative widths of the uncertainty intervals are calculated by dividing the width of the uncertainty interval by the true y and averaging across the L simulations.
In each plot, the x-axis is the metric value and the y-axis denotes different models. The panels represent the different within-cluster sampling schemes. The top two plots are for the fixed-percentage schemes, where n j = ρN j for ρ = 0.1 and ρ = 0.5, j = 1, . . . , J s . The bottom two plots represent the self-weighting samples, with n j = 10 and n j = 50, j = 1, . . . , J s . The colors of the circles represent different first-stage sample sizes J s , J s ∈ {10, 50}.
We now describe the results for each of these three comparisons for the four combinations of outcome type (continuous and binary) and population cluster size model (Poisson and multinomial distributions) as explained in the previous section.
Bayesian methods generally yield more efficient inference than classical estimators, particularly with small number of clusters. For continuous y, the Bayesian models outperform the design-based estimators, both for the Poisson and the multinomially distributed population cluster sizes in Figures 1 and 2, respectively. The differences are rather small when J s = 50 but pronounced when J s = 10. The Hájek estimator has large bias, particularly when the sample is self-weighting, but including auxiliary information as the GREG estimator does greatly reduces the bias. Still, the classical estimators yield unstable results, evident in the high RRMSEs. The Bayesian estimators are preferable here with lower bias and RRMSE, and yield short uncertainty intervals whose coverage rates are close to or above the nominal level. For binary y, there is little difference between the Bayesian methods and the Hájek estimator when the number of sampled clusters is large, J s = 50. This holds for both the Poisson-distributed cluster sizes in Figure 3 and the multinomially distributed cluster sizes in 4. When the number of sampled clusters is small, the Hájek estimator and the Bayesian methods has comparable bias and RRMSE. However, the coverage rates for the Hájek estimator are often below the nominal level, particularly when the sample is not self-weighting (top row of plots).
Both the parametric and nonparametric approaches are statistically valid and have competitive performances. For continuous y the parametric models negbin and lognormal perform comparably to the nonparametric bb with unbiased estimates and similar RRM-SEs in Figures 1 and 2 particularly under large J s , while coverage is generally higher for the parametric models in Figure 2. For binary y, with Poission-distributed cluster sizes the parametric models have a bit higher bias in Figure 3, ranging around 1-1.5%, while for multinomially distributed cluster sizes in Figure 4 the parametric models are less biased than the nonparametric one, especially when the sample is not self-weighting and the number of clusters is small. Coverage rates vary, but are most consistently around or above the nominal level both for the parametric and nonparametric methods. For both continuous and binary y, there is little difference in uncertainty interval lengths between the parametric and nonparametric methods.
Incorrectly omitting cluster sizes as cluster-level predictors-that is, using cluster inds instead of knowsizes-has little impact when y is continuous for either the Poisson or the multinomially distributed population cluster sizes. The bias, RRMSE, and coverage rates for the two methods are similar in both Figures 1 and 2. The differences between cluster inds and knowsizes are minor for binary y as well; cluster inds does not perform appreciably worse than knowsizes in either Figure 3 or 4, the Poisson or the multinomially distributed population cluster sizes. If y and N j are unrelated, it is not necessary to include N j in the model, even under PPS sampling; allowing the regression parameters to vary by cluster as in cluster inds is sufficient for valid inference. In the application study of Section 4, we find that including the cluster sizes as cluster-level predictors will substantially reduce bias and RRMSE with continuous outcome but had subtle difference under binary outcome comparing to the approach only including cluster indicators as random effects models. It's pivotal to account for the two-stage structure comparing to the PPS design. This shows when the sampling design is complex, including two-stage sampling, cluster sampling, PPS and SRS, some design feature could play a bigger role than others. We recommend controlling for all the design features if possible.

Fragile Families Study Application
To evaluate the performance of our method in a more realistic survey context, we use a modified version of the Fragile Families study design in conjunction with a presumed outcome model to implement the finite population inference. We would like to use the Fragile Families sampling frame to illustrate the benefits of Bayesian models accounting for the design features. For convenience, we use the outcome estimation model that is the same as the generation model, which assumption can be released as future extensions.
The Fragile Families study (Reichmann et al., 2001) divided the 77 U.S. cities with 1994 populations of 200,000 or greater into nine strata based on their policy environments and labor markets. Eight of the strata were for cities with extreme values in at least one of the three policy dimensions under consideration (labor markets, child support enforcement, and welfare generosity), and the ninth stratum was for cities that had no extreme values. One city was selected with PPS in each of the eight extreme strata, with a target sample size of 325 births in each city. In the last stratum, eight cities were selected via PPS, with a target sample size of 100 births in each. There was an intermediate stage of selecting hospitals, which we ignore for the paper illustration. We use the Fragile Families study's city population of 77 cities in 1994 as the sampling frame and implement two-stage cluster sampling under PPS.
As a simulation, we use the city population (divided by 100 for computational convenience) as both the measure of size M j and the number of units in the cluster N j , though the ultimate unit of sampling in the study was births and number of births in cities should be accounted for. We exclude the three cities that would be selected with probability one for a total of J = 74 cities. For each unit in the population, we generate an outcome y according to our model in (13). While the original Fragile Families sampling design involves nine strata, we combine them into a single stratum. As in the actual study design, we sample 16 cities with probability proportional to the city population. In each sampled city, we sample either 325 or 100 births, depending on whether the city is a large-or small-sample city, as designated in Reichmann et al. (2001), which results in a self-weighting sample. Figures 5 and 6 show the outputs for when the outcome is continuous and binary, respectively, in terms of relative bias, RRMSE, coverage rates and relative widths of 50% and Hájek: the Hájek estimator; cluster inds: the model with random cluster effects but without the cluster size predictor; knowsizes: the model with known population cluster sizes.
95% uncertainty intervals. The main findings are consistent with the simulation studies. For continuous y in Figure 5, the Bayesian methods (with the exception of negbin) outperform the design-based estimators in terms of RRMSE and uncertainty interval width and are comparable on bias and coverage. The Bayesian methods yield uncertainty intervals that are less than half the width of those based on the design-based methods, with coverage rate that is close to the nominal level. Among the three Bayesian methods, bb and lognormal perform similarly, and both are better than the negbin assumption. The negative binomial population distribution performs poorly with large bias and RRMSE but low coverage rate. Excluding cluster sizes leads to worse performance, with higher bias and RRMSE and longer uncertainty intervals for cluster inds compared to knowsizes.
When y is binary as in Figure 6, we again see that the Bayesian methods yield better results in terms of bias, RRMSE, and coverage than the classical Hájek estimator. The uncertainty intervals of the Hájek estimator are the shortest but are close to those from the Bayesian methods. Comparing the parametric and nonparametric models, bias and RRMSE are lower with higher coverage rates for lognormal than for bb, though coverage rates for both are above the nominal levels. However, lognormal is conservative since the confidence intervals are longer than those for bb. The negbin has the largest bias but shortest uncertainty intervals with comparable RRMSE and coverage among the three approaches. The effects of excluding the cluster sizes are small, with cluster inds having only slightly larger bias and RRMSE than knowsizes.
To further investigate the population distribution of cluster sizes, Figure 7 shows the density plots for 100 cluster sizes drawn from the assumed Poisson and multinomial distributions and the 74 (non-certainty for selection) Fragile Families city populations. From the plots, both the Poisson distribution and the multinomial/Gamma distribution used in the simulation study are different from the population distribution of cluster sizes in the Fragile Family study. The cluster sizes in the Fragile Family study are highly skewed. Hence, in the application, the negative binomial size-biased distribution assumption is not appropriate to depict the cluster size population with poor performance. The Bayesian bootstrap method avoids the parametric assumption and yields robust inference, and the lognormal distribution as the size-biased choice is able to capture the skewness and performs competitively, as shown in Figure 5 and Figure 6. We can modify the parametric assumptions and improve the inference with suitable prior knowledge. negbin: negative binomial size-biased distribution; lognormal: lognormal size-biased distribution; bb: Bayesian bootstrap; Hájek: the Hájek estimator; greg: generalized regression estimator; cluster inds: the model with random cluster effects but without the cluster size predictor; knowsizes: the model with known population cluster sizes. weighted by the odds of that cluster not being selected. While design-based approaches are common in survey inference, variance estimation is often challenging. Current estimation approaches include theoretical approximations  and resampling methods (e.g., Wolter, 2007). In contrast, our integrated approach yields the posterior distribution for the quantities of interest about the finite population, from which variances, uncertainty intervals, and any other functions can easily be computed. The proposal accounts for the design features in modeling and yields designconsistent inference.
The Bayesian methods generally outperform the design-based estimators and improve inference stability, particularly when the number of sampled clusters is small. The performance of the parametric methods negbin and lognormal is comparable to that of the nonparametric Bayesian bootstrap. When extra information about the population cluster sizes is available, for example, from previous years or similar groups, we can incorporate through the informative prior information. Moreover, the parametric methods are straightforward to implement in Stan, which makes them accessible to researchers whose expertise is in areas outside of statistics or programming. The results for parametric and nonparametric methods are more similar when J s = 50 than when J s = 10 in many of the scenarios our simulation study considered. The parametric method is subject to model misspecification especially under small sample. We recommend using the parametric methods as an initial step and perform model diagnostics. An important diagnostic measure is to check whether the population cluster sizes are highly skewed, as in the case of the Fragile Families setup shown in Figure 7. Thus, reasonable prior knowledge of the population distribution of cluster sizes should guide the model choice of parametric or nonparametric approach.
In our study, under binary y, the Bayesian methods were less clearly superior to classical methods in estimating the finite population proportion. One possible reason is that few auxiliary or predictive variables are included in the model. However, when the cluster sizes are highly skewed, as in the Fragile Families case, Bayesian methods perform significantly better, in terms of lower bias and more reasonable coverage, than the classical estimators.
There are several interesting directions in which the current research could be extended. First, our simulation has not considered the case where M j = N j in depth. The natural next step would be to extend the Fragile Families simulation to include the case where the measure of size M j is the city population, but the cluster size N j itself is the total number of births in the city. In doing so, we must make some additional assumptions. So far, we have assumed that we know M j only for the sampled clusters, but what about N j ? If both M j and N j are only available for sampled clusters, we shall predict both M j and N j for the entire population. One idea is to assume that N j is a function of M j and use regression models to predict N j given M j , perhaps the on the log scale to avoid predicting negative cluster sizes and difficulties with cluster sizes ranging over several orders of magnitude. In the Fragile Families study, the correlation between the log of city population M j and log of total births N j is 0.78, so this seems like a promising strategy.
Second, the outcome model can be extended with flexible modeling strategies. To focus on evaluation of different approaches to accounting for the design effect and predicting the nonsampled cluster sizes, for the outcome model, the estimation model we use is the same as the data generation model. In practice, we recommend outcome modeling that is robust against misspecification. Flexible models in the literature can be explored, such as heteroscedasticity assumption, penalized spline regression models, and nonparametric Bayesian models. The multilevel models stabilize estimation via smoothing across clusters. The partial pooling effect can be strengthened with generalized covariance structure, e.g., covariance kernel functions in Gaussian process regression models.
Another direction would be to consider a stratified PPS design as in the original Fragile Families study design. This extension introduces another challenge in that we would need to adjust for the strata structure in our model. For the parametric cluster size models, we would need to partially pool the size parameters (e.g., µ, φ in the negative binomial model, µ, σ in the lognormal) across strata, adding another layer of complexity to the model.
Bayesian approaches are well-equipped to account for the design features in the survey data under complex sampling design through hierarchical modeling. Computational software development, such as the use of Stan, makes modeling approaches enhance the advantage. More methodology developments are necessary to incorporate additional information about the sampling into modeling, such as known population size, paradata and auxiliary variables.