Statistical estimates from survey samples have traditionally been obtained via design-based estimators (Lohr, 2010). In many cases, these estimators tend to work well for quantities such as population totals or means, but can fall short as sample sizes become small. In today’s “information age,” there is a strong demand for more granular estimates. The Small Area Income and Poverty Estimates program (SAIPE) and the Small Area Health Insurance Estimates program (SAHIE) are two examples that rely on American Community Survey (ACS) data, where granularity is essential (Luery, 2011; Bauder ., 2018). Both of these programs require annual estimates at the county level for the entire United States. Many counties exhibit extremely small sample sizes, or even a sample size of zero. In these cases, design-based estimation is inadequate and model-based estimation becomes necessary.
Models for survey data can be either at the area level or the unit level. Area-level models typically use design-based estimators as the response and tend to smooth the estimates in some fashion. These models often use area-level random effects to induce smoothing, and thus a common application is small area estimation (see, for example, Porter . (2015) and the references therein). Rao Molina (2015) provide a recent overview of many of the current area-level models that are available. One issue with area-level models is that estimates at a finer geographic scale may not aggregate to estimates at coarser spatial resolutions, thereby producing inconsistencies.
Unit-level models include individual response values from the survey units as response variables rather than the area-level design-based estimators. The basic unit-level model was introduced byBattese . (1988) in order to estimate small area means. One advantage of unit-level modeling is that the response value can be predicted for all units not contained in the sample, and thus estimates for finite population quantities aggregate naturally. In addition, unit-level models have the potential to yield more precise estimates than area-level models (Hidiroglou You, 2016)
. When modeling survey data at the unit level, the response is often dependent on the sample selection probabilities. This scenario is termedinformative sampling, and it is critical to incorporate the design information into the model in order to avoid biased estimates (Pfeffermann Sverchkov, 2007). Various approaches exist for incorporating an informative design into a model formulation. Little (2012) suggests the use of design variables in the model. For simple survey designs, this may work well, but can become infeasible for complex survey designs. Si . (2015) and Vandendijck . (2016) both use nonparametric regression techniques on the survey weights. These types of techniques do not require any knowledge of the survey design, though they can be difficult to implement in the presence of covariates. Finally, the often used pseudo-likelihood approach (Skinner, 1989; Binder, 1983) exponentially weights each unit’s likelihood contribution by the corresponding survey weight. In this way, the sample data model is adjusted to better match the population distribution. See Parker . (2019) for a recent review.
In addition to the issues that arise due to informative sampling, many variables found within survey data are non-Gaussian in nature, which may induce modeling difficulties. Two examples present in American Community Survey (ACS) data are a binary indicator of health insurance coverage and a count of the number of bedrooms within a household. The SAE setting is often aided by the use of area-level and/or unit-level random effects, which is commonly done using Bayesian hierarchical modeling with a latent Gaussian process (LGP) (Cressie Wikle, 2011; Gelfand Schliep, 2016). In the presence of non-Gaussian data, LGP models lead to non-conjugate full conditional distributions that can be difficult to sample from. Bradley . (2017) provide a solution to this problem by appealing to a class of multivariate distributions that are conjugate with members of the natural exponential family.
We introduce a modeling framework for dealing with unit-level count data under informative sampling by using Bayesian hierarchical modeling to account for complex dependence structures (e.g., in space and time), and relying on the distribution theory provided by Bradley . (2017) for computationally efficient sampling of the posterior distribution. To account for informative sampling, we use a Bayesian pseudo-likelihood (Savitsky Toth, 2016). In Section 2 we introduce and discuss our modeling approach. Section 3 considers a simulation study comparing our methodology to that of two competing estimators. Finally, we provide discussion in Section 4.
2.1 Informative Sampling
Parker . (2019) review current approaches to unit-level modeling under informative sampling. Some of the general approaches include incorporating the design variables into the model (Little, 2012), regression on the survey weights (Si ., 2015; Vandendijck ., 2016), and joint modeling of the response and weights (Pfeffermann Sverchkov, 2007; Novelo Savitsky, 2017). Another general approach is to use a weighted pseudo-likelihood.
Let be an enumeration of the units in the population of interest, and let be the observed, sampled units, selected with probabilities . Let be a variable of interest associated with unit . Our goal is inference on the finite population mean . Suppose a model,
, conditional on a vector of unknown parameters,, holds for the units for . If the survey design is informative, so that the selection probabilities, , are correlated with the response variables, , the model for the nonsampled units will be different from the model for the sampled units, making inference for the finite population mean challenging. Often, the reported survey weights, are used to account for the survey design.
The pseudo-likelihood (PL) approach, introduced by Skinner (1989) and Binder (1983), uses the survey weights to re-weight the likelihood contribution of sampled units. The pseudo-likelihood is given by
where is the response value for unit in the sample . In (1), the vector of model parameters is denoted by and the survey weight for unit is denoted by . For frequentist estimation, the PL can be maximized via maximum likelihood techniques, whereas Savitsky Toth (2016) use the PL in a Bayesian setting for general models. Modeling under a Bayesian pseudo-likelihood induces a pseudo-posterior distribution
represents the weights after being scaled to sum to the sample size. This scaling is done in order to keep the asymptotic amount of information the same as the regular likelihood case, and prevent under-estimation of the standard errors, since the weights act as frequency weights(Savitsky Toth, 2016). It was shown by Savitsky Toth (2016), that the pseudo-posterior distribution converges to the population posterior distribution, justifying the use of the pseudo-posterior distribution for inference on the nonsampled units.
The PL approach is geared towards parameter estimates and not necessarily estimates of finite population quantities. Nevertheless, in this setting (and others), poststratification is a general technique that can be used to create finite population quantity estimates. The general idea is to use a model to predict the response value for all unsampled units in the population, effectively generating a population that can be used to extract any desired estimates. Little (1993) gives an overview of poststratification, whereas Gelman Little (1997) and Park . (2006) develop the idea of poststratification under Bayesian hierarchical models.
2.2 Modeling Non-Gaussian Data
Many of the variables collected from complex surveys are non-Gaussian. For example, binary indicators and count data are both very common in survey data, but cannot be modeled at the unit level under a Gaussian response framework. As such, this can lead to computational issues when dependence structures are introduced.
Bayesian hierarchical modeling is commonly used to model complex dependence structures such as those found in sample surveys. These models often consist of a data stage that models the response, a process stage, and a prior distribution over model parameters. Traditionally, a latent Gaussian process is used to model the process stage; for example see Cressie Wikle (2011). Gelfand Schliep (2016) review the use of Gaussian process modeling in spatial statistics, and Bradley . (2015) develop a general LGP framework that handles multivariate responses as well as complex spatio-temporal dependence structures.
In a Bayesian setting, when the response variable is also Gaussian, Gibbs sampling can be implemented to efficiently sample from the posterior distribution. Unfortunately, when dealing with non-Gaussian data, a Metropolis-Hastings type step may be necessary within the Markov chain Monte Carlo algorithm. Consequently, this algorithm must be tuned and can lead to poor mixing, especially in high-dimensions. Because many survey variables are inherently non-Gaussian, the Gaussian process framework is not ideal in many survey data scenarios.
incorporate new distribution theory to create a set of Bayesian hierarchical models that maintain conjugacy for any response variable contained in the natural exponential family. This includes Poisson, Bernoulli, Binomial, and Gamma random variables, among others, and thus offers a very general modeling framework that maintains computational efficiency.
In this work we consider the distribution theory for Poisson responses specifically, in order to model count survey data. Bradley . (2017) further consider the Negative Binomial case, but state that modeling can be more challenging in this scenario. They suggest that Negative Binomial data may be alternatively modeled as Poisson, and inclusion of random effects can help to model overdispersion.
Each natural exponential family response type is shown to be conjugate with a class of distributions referred to as the conjugate multivariate (CM) distribution. For a Poisson response, the CM distribution is the multivariate log-Gamma (MLG) distribution with probability density function (PDF)
denoted by . The MLG distribution is easy to simulate from using the following steps:
Generate a vector as independent Gamma random variables with shape and rate , for
Bayesian inference with Poisson data and MLG prior distribution also requires simulation from the conditional multivariate log-Gamma distribution (cMLG). Letting , Bradley, Holan Wikle (2018) show that can be partitioned into , where is -dimensional and is -dimensional. The matrix is also partitioned into , where is an matrix and is an matrix. Then
where , and is a normalizing constant. It is also easy to sample from the cMLG distribution when doing Bayesian analysis by using a collapsed Gibbs sampler (Liu, 1994). Bradley . (2017) show that this can be done by drawing , where is sampled from .
2.3 Pseudo-likelihood Poisson Multivariate log-Gamma Model
In order to use the conjugate multivariate distribution theory of Bradley . (2017) in a survey setting for count data under informative sampling, we replace the Poisson likelihood with a survey weighted pseudo-likelihood. Under the unweighted setting, the likelihood contribution to the posterior is proportional to
with representing a vector of response variables, and representing a parameter vector, which will later be modeled using the MLG distribution. The parameter for the Poisson case. This expression is proportional to the product of Poisson densities with natural parameters , . By exponentiating the Poisson likelihood by a vector of weights, , the pseudo-likelihood contribution to the posterior is then proportional to
with representing a Hadamard product, or element-wise multiplication. This is the same form as , where and
, and thus the MLG class of distributions is conjugate with pseudo-likelihoods built upon the Poisson distribution. This is important, as it allows us to use Gibbs sampling with conjugate full conditional distributions in order to sample from the posterior distributions.
Furthermore, Bradley, Holan Wikle (2018) show that the
converges in distribution to a multivariate normal distribution with meanand covariance matrix as the value of approaches infinity. This is convenient as it allows one to effectively use a latent Gaussian process model structure, while still maintaining the computationally benefits of conjugacy offered by the conjugate multivariate distribution theory. Herein, for illustration purposes, we use this type of prior distribution to approximate a latent Gaussian process. However, if desired, one could further model the shape and scale parameters from the MLG prior distribution, which can result in a more flexible shape to the posterior distribution.
We now consider the pseudo-likelihood Poisson multivariate log-Gamma model (PL-PMLG),
where is the th response variable for unit in the sample. This model uses a pseudo-likelihood to account for informative sampling, and is built upon a Poisson response type in order to handle count valued survey data. In this work, the vector corresponds to an incidence vector for which areal unit resides in. As such, the vector acts as area level random effects, which are shared across response types in order to induce multivariate dependence. We note that this model is written for multivariate responses, but we focus only on a univariate example in this work. The parameters act as unit level random effects, and can account for fine scale variation due to missing unit level covariates. Finally, corresponds to fixed effects, for which covariates may or may not be shared across response types. We place log-Gamma priors truncated below at zero (denoted ) on the parameters and . This is done to maintain conjugate full conditional distributions, although other prior distributions could be used here with minimal tuning required as these are low-dimensional parameters deep in the model hierarchy. We set in order to approximate Gaussian prior distributions. We also set in order to create vague prior distributions. However, if prior knowledge on these parameters exists, these values could be adjusted accordingly. The full conditional distributions used for Gibbs sampling can be found in Appendix A.
2.4 Boundary Correction
One technical issue that arises when using a conjugate multivariate hierarchical modeling framework concerns data that are observed on the boundary of their support (i.e. zero counts for Poisson data). When zero counts are observed with Poisson data, the result is a full conditional distribution with a shape parameter of zero which is not well defined. Because the conjugate multivariate framework was only recently developed, there is relatively little literature on handling these boundary issues; however Bradley . (2017) suggest using adjusted data, , by adding a small constant. This can work in many cases, depending on the dataset and the value of , but is effectively sampling from an approximation to the posterior distribution.
Rather than sample from an approximate distribution, we use importance sampling to sample from the true posterior distribution, similar to the work of Kim . (1998). In this case, the importance weights are proportional to the ratio of the adjusted pseudo-likelihood to the true pseudo-likelihood. However, with large sample sizes, the adjusted pseudo-likelihood can diverge from the true pseudo-likelihood. To this effect, we run a pilot chain (using 100 iterations) to find the average ratio of the true log-pseudo-likelihood to the adjusted log-pseudo-likelihood. We then scale the weights in the adjusted pseudo-likelihood by this average ratio. This has the effect of centering the adjusted pseudo-likelihood around the true pseudo-likelihood. The importance weights, taken at each iteration of the Gibbs sampler, are then proportional to
where represents the scaled survey weight after multiplying by the average ratio mentioned above. We found that for the constant a value of one or two was ideal, as it minimized the extent of the divergence from the true pseudo-likelihood to the approximate one.
3 Empirical Simulation Study
The American Community Survey (ACS) is an ongoing survey, with approximately 3.5 million households sampled annually, that is critical for informing how federal funds should be allocated. Although the complete microdata is not available to the public, public use microdata samples (PUMS) are available. PUMS only contain geographic indicators at the public use microdata area level (PUMA), which are aggregated areas such that each contains at least a population of 100,000 people. For this survey as well as others, vacant houses can pose a challenge when conducting the survey. Bradley . (2017) use a simple random sample of ACS PUMS data within a single PUMA to predict housing vacancies by modeling the number of people per household as a Poisson random variable. This work illustrates the capacity of unit-level models to predict housing vacancies, however, because they used simple random sampling within a single PUMA, the methodology cannot be applied in an informative sampling context for SAE.
We construct an empirical simulation study to illustrate how the PL-PMLG can be used to create small area estimates of the number of housing vacancies. Using the state of Alabama, we treat the entire 2017 PUMS housing dataset as our population (or “truth”). This dataset contains roughly 22,500 observations across 34 different PUMAs. We further subsample this data using the Midzuno probability proportional to size method (Midzuno, 1951) within the ‘sampling’ R package (Tillé Matei, 2016) , which we use to create our estimates. We then compare these estimates to the truth.
In addition to comparing the PL-PMLG to a direct estimator, we also wish to compare to another model based estimator. In this scenario, many of the direct estimators are equal to zero, which makes area-level modeling prohibitively difficult. Instead, because count data are often modeled as Gaussian on the log scale, we compare to a unit level model taking this approach. Because the data contains zero counts, a small constant, , must be added to the data before taking the log transformation, and this transformation is undone when predictions are made. The full model hierarchy, which we call the Gaussian Approximation model (GA), is
where we use the vague prior distribution , and . We again use a pseudo-likelihood approach here in order to account for informative sampling. The rest of the model consists of fairly standard Bayesian mixed effects regression. We tested the value fixed over the values of (0.1, 1, 5), and found that yielded substantially lower MSE and bias for this example, which is what we present here.
For this simulation, we take a sample size of 5,000 from the PUMS data with probability proportional to
(i.e., probability inversely proportional to the original probability of selection). We show that sampling this way induces informativeness by comparing to the unweighted version of our model. Our fixed effects consist of an intercept, and the number of bedrooms in the household, which we treat as a categorical variable. We calculate the Horvitz-Thompson estimate (direct estimate) as well as the two model-based estimates. Finally, we repeat the process 50 times in order to compare MSE and absolute bias. For the PL-PMLG, unweighted PMLG (UW-PMLG) and GA estimators, we used Gibbs sampling for 2,000 iterations, discarding the first 1,000 as burn-in. Convergence was assessed visually through traceplots of the sample chains, and no lack of convergence was detected. We also compare to a Horvitz-Thompson direct estimator, with Hajek variance estimates using themase package in R (McConville ., 2018).
A summary of the simulation results can be found in Table 1, where we compare the MSE and absolute bias of the PUMA level estimates for the total number of vacant housing units. The GA model does not provide a reduction in MSE compared to the direct estimator; however, the unweighted and PL-PMLG models do (12% and 49% respectively). Additionally, the absolute bias for the PL-PMLG is substantially lower than the GA and unweighted models. The significant reduction in MSE and bias comparing the PL-PMLG and UW-PMLG models indicates that there was an informative design, and the PL approach helps to account for this design. We also show the point estimates from a randomly chosen single run of the simulation under each estimator in Figure 1. All of the estimators seem to capture the same general spatial trend, however the PL-PMLG estimator seems to most closely resemble the truth. As a final comparison, we plot the standard error of the estimates averaged across the 50 simulations on the log scale in Figure 2
. To construct this figure, we compute a standard error of the estimate under each approach, for each of the 50 simulated datasets. For the model-based estimates, this standard error is the posterior predictive standard deviation. We then average these standard errors across the simulated data sets, in order to illustrate the expected uncertainty associated with each reported estimate. In some cases, the standard error of the direct estimate could not be obtained due to a point estimate of zero, in which case they have been removed from the average. As expected, the standard errors are dramatically lower for the model-based estimators than the direct estimator. In general the GA standard errors are slightly lower than the PL-PMLG, however the GA exhibits much higher MSE due to the increased bias, as evidenced by Table1. Thus, the PL-PMLG appears to be a superior estimator overall.
There is a strong need for unit-level models that can handle survey data. Accounting for informative sampling design and modeling non-Gaussian data types are two of the biggest challenges in this setting. In this work, we present a new method for modeling count data while accounting for informative sampling. This method can be used for SAE as well as for more general modeling purposes. Our method relies on conjugate multivariate distribution theory, and we show that conjugacy is maintained when using a psuedo-likelihood approach to account for the survey design. We also extend the work of Bradley . (2017) to handle the issue of zero counts through importance sampling.
Our approach is illustrated on a simulation study built upon public-use ACS data. This is a count data example where area-level models are not feasible and Gaussian models are not appropriate. Furthermore, this is an example where direct estimators are not useful due to excessively large MSE and standard errors. Our PL-PMLG approach is able to accurately estimate population quantities based on count variables while still maintaining computational efficiency.
There still remains further work to be done in the area of non-Gaussian survey data. Other data types such as binary random variables are prevalent and should be considered (Bauder ., 2018; Luery, 2011). The conjugate multivariate framework offered by Bradley . (2017) has the potential to fit these types of data, although the boundary value issue may pose a computational challenges. Our importance sampling approach works well in the Poisson case, but a more general solution may be attainable. Finally, non-Gaussian data should be explored in regards to other solutions to the informative sampling problem. The pseudo-likelihood approach may be one of the most popular approaches to informative sampling, but other methods exist and may yield additional gains in terms of precision.
This report is released to inform interested parties of ongoing research and to encourage discussion of work in progress. The views expressed are those of the authors, and not those of the U. S. Census Bureau. This research was partially supported by the U.S. National Science Foundation (NSF) and the U.S. Census Bureau under NSF grant SES-1132031, funded through the NSF-Census Research Network (NCRN) program and NSF SES-1853096. Partial support of this work through the U.S. Census Bureau Dissertation Fellowship Program is gratefully acknowledged.
Appendix A Full Conditional Distributions for the PL-PMLG Model
a.1 Random Effects
a.2 Fixed Effects
a.3 Variance Parameters
- Battese . (1988) bat88Battese, GE., Harter, RM. Fuller, WA. 1988. An error-components model for prediction of county crop areas using survey and satellite data An error-components model for prediction of county crop areas using survey and satellite data. Journal of the American Statistical Association8340128–36.
- Bauder . (2018) bau18Bauder, M., Luery, D. Szelepka, S. 2018. Small area estimation of health insurance coverage in 2010 – 2016 Small area estimation of health insurance coverage in 2010 – 2016 . Small Area Methods Branch, Social, Economic, and Housing Statistics Division, U. S. Census Bureau. https://www2.census.gov/programs-surveys/sahie/technical-documentation/methodology/2008-2016-methods/sahie-tech-2010-to-2016.pdf
- Binder (1983) bin83Binder, DA. 1983. On the variances of asymptotically normal estimators from complex surveys On the variances of asymptotically normal estimators from complex surveys. International Statistical Review513279–292.
- Bradley . (2015) brad15Bradley, JR., Holan, SH. Wikle, CK. 2015. Multivariate spatio-temporal models for high-dimensional areal data with application to Longitudinal Employer-Household Dynamics Multivariate spatio-temporal models for high-dimensional areal data with application to longitudinal employer-household dynamics. The Annals of Applied Statistics941761–1791.
- Bradley . (2017) brad17Bradley, JR., Holan, SH. Wikle, CK. 2017. Bayesian Hierarchical Models with Conjugate Full-Conditional Distributions for Dependent Data from the Natural Exponential Family Bayesian hierarchical models with conjugate full-conditional distributions for dependent data from the natural exponential family. arXiv preprint arXiv:1701.07506.
- Bradley, Holan Wikle (2018) brad18Bradley, JR., Holan, SH. Wikle, CK. 2018. Computationally Efficient Multivariate Spatio-Temporal Models for High-Dimensional Count-Valued Data (with Discussion) Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data (with discussion). Bayesian Analysis131253–310.
- Bradley, Wikle Holan (2018) brad18bBradley, JR., Wikle, CK. Holan, SH. 2018. Spatio-Temporal Models for Big Multinomial Data using the Conditional Multivariate Logit-Beta Distribution Spatio-temporal models for big multinomial data using the conditional multivariate logit-beta distribution. arXiv preprint arXiv:1812.03555.
- Cressie Wikle (2011) cress11Cressie, N. Wikle, CK. 2011. Statistics for Spatio-Temporal Data Statistics for spatio-temporal data. John Wiley & Sons.
- Gelfand Schliep (2016) gelf16Gelfand, AE. Schliep, EM. 2016. Spatial statistics and Gaussian processes: A beautiful marriage Spatial statistics and Gaussian processes: A beautiful marriage. Spatial Statistics1886–104.
Gelman Little (1997)
gelman97Gelman, A. Little, TC.
Poststratification into many categories using hierarchical logistic regression Poststratification into many categories using hierarchical logistic regression.Survey Methodology23127–135.
- Hidiroglou You (2016) hid16Hidiroglou, MA. You, Y. 2016. Comparison of unit level and area level small area estimators Comparison of unit level and area level small area estimators. Survey Methodology4241–61.
- Kim . (1998) kim98Kim, S., Shephard, N. Chib, S. 1998. Stochastic volatility: likelihood inference and comparison with ARCH models Stochastic volatility: likelihood inference and comparison with arch models. The review of economic studies653361–393.
- Little (1993) little93Little, RJ. 1993. Post-stratification: a modeler’s perspective Post-stratification: a modeler’s perspective. Journal of the American Statistical Association884231001–1012.
- Little (2012) little12Little, RJ. 2012. Calibrated Bayes, an alternative inferential paradigm for official statistics Calibrated bayes, an alternative inferential paradigm for official statistics. Journal of Official Statistics283309.
- Liu (1994) liu94Liu, JS. 1994. The collapsed Gibbs sampler in Bayesian computations with applications to a gene regulation problem The collapsed gibbs sampler in bayesian computations with applications to a gene regulation problem. Journal of the American Statistical Association89427958–966.
- Lohr (2010) lohr09Lohr, S. 2010. Sampling: Design and Analysis, 2nd ed Sampling: Design and analysis, 2nd ed. Boston: Brooks/Cole.
- Luery (2011) luery11Luery, DM. 2011. Small area income and poverty estimates program Small area income and poverty estimates program. Proceedings of 27th SCORUS Conference Proceedings of 27th scorus conference ( 93–107).
- McConville . (2018) maseMcConville, K., Tang, B., Zhu, G., Cheung, S. Li, S. 2018. mase: Model-Assisted Survey Estimation mase: Model-assisted survey estimation . https://cran.r-project.org/package=mase
- Midzuno (1951) midz51Midzuno, H. 1951. On the sampling system with probability proportionate to sum of sizes On the sampling system with probability proportionate to sum of sizes. Annals of the Institute of Statistical Mathematics3199–107.
- Novelo Savitsky (2017) novelo17Novelo, LL. Savitsky, T. 2017. Fully Bayesian Estimation Under Informative Sampling Fully Bayesian estimation under informative sampling. arXiv preprint arXiv:1710.00019.
- Park . (2006) park06Park, DK., Gelman, A. Bafumi, J. 2006. State-level opinions from national surveys: Poststratification using multilevel logistic regression State-level opinions from national surveys: Poststratification using multilevel logistic regression. Public Opinion in State Politics.
- Parker . (2019) par19Parker, PA., Janicki, R. Holan, SH. 2019. Unit Level Modeling of Survey Data for Small Area Estimation Under Informative Sampling: A Comprehensive Overview with Extensions Unit level modeling of survey data for small area estimation under informative sampling: A comprehensive overview with extensions. arXiv preprint arXiv:1908.10488.
- Pfeffermann Sverchkov (2007) pfe07Pfeffermann, D. Sverchkov, M. 2007. Small-area estimation under informative probability sampling of areas and within the selected areas Small-area estimation under informative probability sampling of areas and within the selected areas. Journal of the American Statistical Association1024801427–1439.
- Porter . (2015) porter15Porter, AT., Wikle, CK. Holan, SH. 2015. Small area estimation via multivariate Fay–Herriot models with latent spatial dependence Small area estimation via multivariate Fay–Herriot models with latent spatial dependence. Australian & New Zealand Journal of Statistics57115–29.
- Rao Molina (2015) rao15Rao, JNK. Molina, I. 2015. Small Area Estimation Small area estimation. Hoboken, New JerseyWiley.
- Savitsky Toth (2016) sav16Savitsky, TD. Toth, D. 2016. Bayesian estimation under informative sampling Bayesian estimation under informative sampling. Electronic Journal of Statistics1011677–1708.
- Si . (2015) si15Si, Y., Pillai, NS. Gelman, A. 2015. Bayesian nonparametric weighted sampling inference Bayesian nonparametric weighted sampling inference. Bayesian Analysis103605–625.
Domain means, regression and multivariate analysis Domain means, regression and multivariate analysis.CJ. Skinner, D. Holt TMF. Smith (), Analysis of Complex Surveys Analysis of complex surveys ( 80 - 84). ChichesterWiley.
- Tillé Matei (2016) sampTillé, Y. Matei, A. 2016. sampling: Survey Sampling sampling: Survey sampling . https://CRAN.R-project.org/package=sampling R package version 2.8
- Vandendijck . (2016) van16Vandendijck, Y., Faes, C., Kirby, RS., Lawson, A. Hens, N. 2016. Model-based inference for small area estimation with sampling weights Model-based inference for small area estimation with sampling weights. Spatial Statistics18455–473.