1 Introduction
Statistical estimates from survey samples have traditionally been obtained via designbased 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, designbased estimation is inadequate and modelbased estimation becomes necessary.
Models for survey data can be either at the area level or the unit level. Arealevel models typically use designbased estimators as the response and tend to smooth the estimates in some fashion. These models often use arealevel 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 arealevel models that are available. One issue with arealevel models is that estimates at a finer geographic scale may not aggregate to estimates at coarser spatial resolutions, thereby producing inconsistencies.
Unitlevel models include individual response values from the survey units as response variables rather than the arealevel designbased estimators. The basic unitlevel model was introduced by
Battese . (1988) in order to estimate small area means. One advantage of unitlevel 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, unitlevel models have the potential to yield more precise estimates than arealevel 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 termed
informative 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 pseudolikelihood 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 nonGaussian 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 arealevel and/or unitlevel 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 nonGaussian data, LGP models lead to nonconjugate 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 unitlevel 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 pseudolikelihood (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 Methodology
2.1 Informative Sampling
Parker . (2019) review current approaches to unitlevel 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 pseudolikelihood.
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 pseudolikelihood (PL) approach, introduced by Skinner (1989) and Binder (1983), uses the survey weights to reweight the likelihood contribution of sampled units. The pseudolikelihood is given by
(1) 
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 pseudolikelihood induces a pseudoposterior distribution
where
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 underestimation of the standard errors, since the weights act as frequency weights
(Savitsky Toth, 2016). It was shown by Savitsky Toth (2016), that the pseudoposterior distribution converges to the population posterior distribution, justifying the use of the pseudoposterior 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 NonGaussian Data
Many of the variables collected from complex surveys are nonGaussian. 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 spatiotemporal 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 nonGaussian data, a MetropolisHastings 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 highdimensions. Because many survey variables are inherently nonGaussian, the Gaussian process framework is not ideal in many survey data scenarios.
Bradley . (2017); Bradley, Holan Wikle (2018); Bradley, Wikle Holan (2018)
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 logGamma (MLG) distribution with probability density function (PDF)
(2) 
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

Let

Let

Then
Bayesian inference with Poisson data and MLG prior distribution also requires simulation from the conditional multivariate logGamma 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
with density
(3) 
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 Pseudolikelihood Poisson Multivariate logGamma 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 pseudolikelihood. 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 pseudolikelihood contribution to the posterior is then proportional to
with representing a Hadamard product, or elementwise multiplication. This is the same form as , where and
, and thus the MLG class of distributions is conjugate with pseudolikelihoods 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 mean
and 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 pseudolikelihood Poisson multivariate logGamma model (PLPMLG),
(4) 
where is the th response variable for unit in the sample. This model uses a pseudolikelihood 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 logGamma 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 lowdimensional 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 pseudolikelihood to the true pseudolikelihood. However, with large sample sizes, the adjusted pseudolikelihood can diverge from the true pseudolikelihood. To this effect, we run a pilot chain (using 100 iterations) to find the average ratio of the true logpseudolikelihood to the adjusted logpseudolikelihood. We then scale the weights in the adjusted pseudolikelihood by this average ratio. This has the effect of centering the adjusted pseudolikelihood around the true pseudolikelihood. 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 pseudolikelihood 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 unitlevel 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 PLPMLG 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 PLPMLG 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 arealevel 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
(5) 
where we use the vague prior distribution , and . We again use a pseudolikelihood 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 HorvitzThompson estimate (direct estimate) as well as the two modelbased estimates. Finally, we repeat the process 50 times in order to compare MSE and absolute bias. For the PLPMLG, unweighted PMLG (UWPMLG) and GA estimators, we used Gibbs sampling for 2,000 iterations, discarding the first 1,000 as burnin. Convergence was assessed visually through traceplots of the sample chains, and no lack of convergence was detected. We also compare to a HorvitzThompson direct estimator, with Hajek variance estimates using the
mase 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 PLPMLG models do (12% and 49% respectively). Additionally, the absolute bias for the PLPMLG is substantially lower than the GA and unweighted models. The significant reduction in MSE and bias comparing the PLPMLG and UWPMLG 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 PLPMLG 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 modelbased 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 modelbased estimators than the direct estimator. In general the GA standard errors are slightly lower than the PLPMLG, however the GA exhibits much higher MSE due to the increased bias, as evidenced by Table
1. Thus, the PLPMLG appears to be a superior estimator overall.Estimator  MSE  Abs. Bias 

Direct  2250  3.5 
GA  2526  33.9 
PLPMLG  1151  23.5 
UWPMLG  1983  32.7 
4 Discussion
There is a strong need for unitlevel models that can handle survey data. Accounting for informative sampling design and modeling nonGaussian 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 psuedolikelihood 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 publicuse ACS data. This is a count data example where arealevel 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 PLPMLG 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 nonGaussian 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, nonGaussian data should be explored in regards to other solutions to the informative sampling problem. The pseudolikelihood 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.
Acknowledgments
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 SES1132031, funded through the NSFCensus Research Network (NCRN) program and NSF SES1853096. Partial support of this work through the U.S. Census Bureau Dissertation Fellowship Program is gratefully acknowledged.
Appendix A Full Conditional Distributions for the PLPMLG Model
a.1 Random Effects
a.2 Fixed Effects
a.3 Variance Parameters
References
 Battese . (1988) bat88Battese, GE., Harter, RM. Fuller, WA. 1988. An errorcomponents model for prediction of county crop areas using survey and satellite data An errorcomponents 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/programssurveys/sahie/technicaldocumentation/methodology/20082016methods/sahietech2010to2016.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 spatiotemporal models for highdimensional areal data with application to Longitudinal EmployerHousehold Dynamics Multivariate spatiotemporal models for highdimensional areal data with application to longitudinal employerhousehold dynamics. The Annals of Applied Statistics941761–1791.
 Bradley . (2017) brad17Bradley, JR., Holan, SH. Wikle, CK. 2017. Bayesian Hierarchical Models with Conjugate FullConditional Distributions for Dependent Data from the Natural Exponential Family Bayesian hierarchical models with conjugate fullconditional 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 SpatioTemporal Models for HighDimensional CountValued Data (with Discussion) Computationally efficient multivariate spatiotemporal models for highdimensional countvalued data (with discussion). Bayesian Analysis131253–310.
 Bradley, Wikle Holan (2018) brad18bBradley, JR., Wikle, CK. Holan, SH. 2018. SpatioTemporal Models for Big Multinomial Data using the Conditional Multivariate LogitBeta Distribution Spatiotemporal models for big multinomial data using the conditional multivariate logitbeta distribution. arXiv preprint arXiv:1812.03555.
 Cressie Wikle (2011) cress11Cressie, N. Wikle, CK. 2011. Statistics for SpatioTemporal Data Statistics for spatiotemporal 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.
1997.
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. Poststratification: a modeler’s perspective Poststratification: 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: ModelAssisted Survey Estimation mase: Modelassisted survey estimation []. https://cran.rproject.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. Statelevel opinions from national surveys: Poststratification using multilevel logistic regression Statelevel 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. Smallarea estimation under informative probability sampling of areas and within the selected areas Smallarea 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.

Skinner (1989)
ski89Skinner, CJ.
1989.
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.Rproject.org/package=sampling R package version 2.8
 Vandendijck . (2016) van16Vandendijck, Y., Faes, C., Kirby, RS., Lawson, A. Hens, N. 2016. Modelbased inference for small area estimation with sampling weights Modelbased inference for small area estimation with sampling weights. Spatial Statistics18455–473.
Comments
There are no comments yet.