{'fraction_non_alphanumeric': 0.0317282985141325, 'fraction_numerical': 0.011060216735560272, 'mean_word_length': 4.512007389162561, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatio-Temporal Mixed Models to Predict Coverage Error Rates at Local Areas

Despite of the great efforts during the censuses, occurrence of some nonsampling errors such as coverage error is inevitable. Coverage error which can be classified into two types of under-count and overcount occurs when there is no unique bijective (one-to-one) mapping between the individuals from the census count and the target population -- individuals who usually reside in the country (de jure residences). There are variety of reasons make the coverage error happens including deficiencies in the census maps, errors in the field operations or disinclination of people for participation in the undercount situation and multiple enumeration of individuals or those who do not belong to the scope of the census in the overcount situation. A routine practice for estimating the net coverage error is subtracting the census count from the estimated true population, which obtained from a dual system (or capture-recapture) technique. Estimated coverage error usually suffers from significant uncertainty of the direct estimate of true population or other errors such as matching error. To rectify the above-mentioned problem and predict a more reliable coverage error rate, we propose a set of spatio-temporal mixed models. In an illustrative study on the 2010 census coverage error rate of the U.S. counties with population more than 100,000, we select the best mixed model for prediction by deviance information criteria (DIC) and conditional predictive ordinate (CPO). Our proposed approach for predicting coverage error rate and its measure of uncertainty is a full Bayesian approach, which leads to a reasonable improvement over the direct coverage error rate in terms of mean squared error (MSE) and confidence interval (CI) as provided by the U.S. Census Bureau.


Introduction
Coverage error can lead to serious problems at both large and small area levels. It might affect policy-making and planning decisions of the country (see Freedman and Navidi (1992)) and produce serious inequity in the administration of federal programs, especially programs which distribute funds based on the statistical formulas rely on the census count. As each year billions of dollars are distributed among areas through such formula grant programs, the concern is understandable (Slater (1990)). Some specific consequences of census undercoverage which might be occured more (c.f. Nations (2010)) is underestimation in postcensal estimates, and overestimation in the use of census data as the denominator in ratio calculations. In addition, it removes part of the target population in surveys where census figures are sampling frames, and it brings bias in the estimators that the census count is the basis for control totals (Slater (1990)).
In order to evaluate the census performance, countries estimate coverage error after the censuses. a negative number indicates the net overcount. By letting T and C to be true population and census count, respectively; we can define T as C + N (Dolson (2010)). Most of countries after the census count estimate the true population (T ) according to one of the existing methods. Some of these methods include demographic analysis, reverse record check methodology, and capture recapture technique. An interested reader might consult Kerr (1998) andU.S. Census Bureau (1985) to get more information about different methods of true population estimation. Then, the net coverage error can be estimated by subtracting census count from the true population estimation, i.e.N =T − C, and eventually net coverage error rate could be defined as (N /T ) × 100.
The capture recapture technique so-called dual system (DS) is one of the routine methods designed for estimating the true population. It is usually based on the information of census and post-enumeration survey (PES)-a probability survey conducted after the census. Its fundamental was firstly set by the United States Census Bureau (USCB) USCB (1985) in the 1950s. In this survey, a representative sample of the census population is enumerated. Homeless persons and institutional households -such as nursing homes, prisons, and dormitories-usually are not counted in the PES as assessing their errors is difficult and selection process is commonly based on housing units/households. Then all people in the selected households or housing units are enumerated.
Estimating the coverage error at the small areas is as important as the large areas such as states or municipalities (Bryce (1990)). The sample size of PES might not be sufficient enough to protect a reliable direct estimate for the coverage error from DS. For instance, if we take a look at the U.S. Census Bureau web page https://www.census.gov/coverage_measurement/post-enumeration_ surveys/2010_results.html, we can realize that the net coverage error estimation even for counties with 100,000+ total population is not reliable enough and has large root mean squared error because of small sample size as one of the crucial reasons. Here it might be of interest to work on some small area modeling to provide more reliable estimates of coverage error rates for the counties. Therefore, the main objective of this document is proposing a new model based estimator for predicting the coverage error rate at the small area levels as there is a growing need of statistics for these areas during the recent years. In our model, we assume the small area is nested in a large area, which has large enough sample size and moderate precise direct estimate of true population and coverage error rate as a consequence. Our method can be used to predict net coverage error rate at the county level which is nested in the state level for the U.S. concern. The proposed model is a new evolution of spatio-temporal mixed model which will be intensely discussed further.
For model building, we borrow strength over three different sources: time-related function of covariates, time smoother, and space smoother. We show later how to make inferential statements that incorporate uncertainty about all of the parameters in the model. That is we do not rely on the simpler empirical Bayes methods that frequently do not incorporate all of the sources of variation, and we rely on a hierarchical Bayesian model. The rest of paper is organized as follows. In Section 2, we exclusively give an overview and explanation through our spatio-temporal mixed models which are well-connected to our salient goals. In Section 3, we give an illustrative study related to the assessment of coverage error rates at the U.S. counties level selected from those with population 100,000+. Finally, we go over the limitations and advantages of models and studies which help us to point the way to where we go later.

Spatio-Temporal Mixed Models
Usually statistical inferences based on the direct estimates under the usual design-based mode of inference do not depend on the validity of the model. However, the design based inference becomes problematic when the sample sizes in subareas are small. In this situation, model based inferences routinely used to produce indirect more reliable estimates. Cressie (1992) proposed a Bayesian model to predict the adjustment factor defined by T /C at the 51 U.S. state levels based on some covariates. Another significant study of Cressie was done in 1990 to predict the adjustment factors based on spatial dependency among states. In his paper, the spatial dependency was based on 700-mile neighbourhoods of states (i.e. : I {d ij ≤700 miles} = 1), where d ij is the distance between the centers of gravity of the i-th and j-th states. The pragmatic role of this spatial model is that an important explanatory variable might be missed from the model (c.f. Cressie (1990)), and a linear relationship of covariates and adjustment factors might be in fact more complicated.
In addition, Cressie limited his work to the space not to the time direction dependency, and he did not go further to the low enough levels such as counties. Finally, he considered the empirical Bayes by providing the naive mean squared error which does not take into account the variability of unknown parameters estimation. On the other hand, the spatial dependency in our model take into account the privilege of combining the cogniguity and distance. Moreover, time could be an important factor in our model and we have the felexibility of nonlinear function of covariates in our model.
More specifically, our objective here is to discuss how to combine the relevant auxiliary informa-tion with time and space in a model to derive a prediction for the coverage error rate at local areas such as counties nested in the states for the U.S. context. Here we believe the sparsity of available data for the small areas-specifically for the next section-require space-smoother and time-smoother (Mercer and Clark (2015)), and we provide a mixed model containing these terms typically called spatio-temporal mixed model. The time-series part of the model looks to see how a particular value is influenced by its past values, and the spatial part looks to see how a particular value is influenced by its neighbouring values.
For the time-series part of the model, a temporal function of some socioeconomic variables such as family status, income, and crime rate more likely related to coverage error rate can be considered. This needs to have an appropriate time-smoother term defined by γ t ∼ IRW (σ 2 γ ) in the model where t and IRW respectively stand for time and intrinsic random-walk notations and σ 2 γ is an unknown scale-parameter term. For the spatial part of the model, we choose the intrinsic conditional auto-regressive (ICAR) structure over other spatial dependency structures as ICAR has primarily better ability to capture relatively smooth spatial dependency. We define the spacesmoother term symbolically by φ s ∼ ICAR(σ 2 u , ρ) which follows φ s |σ 2 u , ρ ∼ N (0, σ 2 u (I − ρW ) −1 M ) where W is a weighting matrix and (σ 2 u , ρ) t is a vector of unknown parameters. The matrix M will be defined later after introducing W as it is a linear map of W .
Our weight matrix W is a function of contiguity δ ij and distance d ij where i and j are similar types of areas while area i-th is the objective one. Contiguity is both rooks and bishops of a county (the areas sharing edge and boundary with the objective county). As all the events occur in R and we need to take into account the nature of under coverage and over coverage in our structure, we propose the following set-up for the weight matrix: For the conventional we assume the two parameters a ≡ b ≡ 1, and we exclude "self influence" by assuming that w ii = 0 for all i = 1, ..., n.
Here, we have considered the Euclidean distance which can be calculated by the formula where ξ i and ζ i are the geographic coordinates (latitude and longitude) of the centroid of each county. More specifically, contiguity can be defined as where ∂ i is a symbol for the neighbourhood of i-th area and y i and y j are coverage error rates for area i and j, respectively, which requires to be in the same direction by y i .y j > 0. The remained matrix M in the ICAR structure is a diagonal matrix M = diag(1/ k j=1 w ij ). In our approach, the matrix (I − ρW ) is invertible and thus fulfills the non-singularity requirement of a spatial auto-regressive equation; we need restrictions on the value of ρ. If W is standardised to have row sums of unity then these restrictions effectively amount to −1 ≤ ρ ≤ 1 (Cliff and Ord (1981)). We might realize that the specific structure that we are interested in cannot be effectively captured from the spatial part of the model; therefore, we can put it in the prior part.
Under some mild conditions, we introduce the following prior for the spatial autocorrelation ρ: This prior was proposed by Kass and Wasserman (1996) which helps us to get more mass at the end-points to maximumly benefit from the spatial correlation. Now, we can put pieces together and introduce our model. As we mentioned earlier, our objective here is predicting the y i (∀i = 1, ..., m) where m is the total number of small areas via a suitable spatio-temporal mixed model. We assume normal distribution for y ts il as follows The hierarchical Bayesian space-time mixed model containing all of terms is y ts il = µ + λ l + β t il + γ t + φ s + δ ts + ts il , ∀i = 1, ..., m l ; l = 1, ..., L.
Here we assume µ (intercept) and λ l (large area effect) are fixed. Also β t il ∼ N (f (x (0,t] il β), τ 2 ) represents a temporal function of vector of covariates based on the pre-specified time-interval (0,t] mildly correlated to y ts il following a normal distribution where τ 2 is unknown. Time smoother and space smoother are γ t ∼ IRW (σ 2 γ ) and φ s ∼ ICAR(σ 2 u , ρ). The space time interaction term is δ ts ∼ iid N (0, σ 2 δ ). Finally, ts il is the sampling design error normally distributed by N (0,V ts il,DS ). The candidate mixed models that we consider in this paper are given in Table 1. A collective vector of unknown parameters is Θ ≡ (µ, λ l , β, τ 2 , σ 2 γ , σ 2 u , σ 2 δ ) t . We assume π(µ, λ l ) ∝ 1 and π(β, logτ 2 ) ∝ c, constant. In addition, we have σ 2 γ , σ 2 u , σ 2 δ ∼ Inv − Gamma(a, b) with a = b = 0.025. After providing the full conditional distributions of unknown parameters, we can generate samples by MCMC numerical integration technique, in particular the Gibbs sampler to obtain the hierarchical Bayes estimate of θ ts il and its corresponding Bayes risk from the full conditional distribution [θ ts il | rest of parameters and data].  model, we conduct a spatial data analysis to model the spatial dependency.
To ensure that we have proposed the right spatial matrix, we discuss techniques for exploring spatial dependency through variogram equivalent to second order properties of spatial dependency.
To do so, an exploratory analysis is to plot (y i − y j ) 2 for all possible permutation pairs (i, j) against the distance between them d ij shown in Figure 2 1 . Figure 2 backs up any informal inferences drawn from the variogram, and we have given separate plots for the undercount and overcount from the overall coverage error to detect any possibilities of variogram strength. The cloud structure of panels in Figure 2 confirms the absence of spatial continuity. Therefore the second order effects do not carry out a pivotal role in the spatial matrix weight. In this case of not such a dependency, we consider the neighbourhood and contiguity weight matrix proposed in the former section. As we can see from Figure 1, some counties are isolated, and they do not have any neighbours. This further problem for which there is no current best advice is how to proceed if some areal entities have no neighbours. When measures of autocorrelation were developed, it was generally assumed that all entities would have neigbours. These areas do not accept spatial weights because of no neighbour entities. So we have decided to remove counties with column sum or row sum zero from the spatial matrix 2 to have a square matrix. This might make the number of neigbours of each county small. Finally, 453 counties which nested in 44 unique states remained in our study shown in Figure 3.
To check the potential interaction between spatial weight matrix and percentage of net coverage error rate for counties, we consider Moran's I statistic which perhaps is the most common global test for our purpose (Bivand and Gomez-Rubio (2008)). As we can see from the Moran's I results (Table   2), the correlation is low since we have small number of neighbours for each county. However, we are currently investigating the possibilities of other spatial weighting matrix than the one proposed in 1 We consider the latitude and longitude of counties from the 2010 Census Gazetteer Files of the Census Bureau https://www.census.gov/geo/maps-data/data/gazetteer2010.html 2 For spatial matrix, we consider the disjunction logic "∨" between columns and rows, meaning that a county could be removed from the spatial matrix iff its sum row or sum column or both are zero in the initial spatial matrix.  Although we include here percentage of net coverage error rates for only 2010 in the spatial part of the model, it is well known that this quantity is of course highly correlated with function of covariates from year to year. Therefore, we have put a five-year time interval 2005 up to 2010 for the covariates in the model to effectively capture the strength of covariates in the i-th location.
The temporal function of covariates contain minority percentage, crime rate, poverty percentage, and percentage with language difficulty coming from QuickFacts UNITED STATES http://www. census.gov/quickfacts/table/PST045215/00.
QuickFacts data are derived from: population estimates, American community survey, census of population and housing, current population survey, small area health insurance estimates, small area income and poverty estimates, state and county housing unit estimates, county business patterns, nonemployer statistics, economic census, and survey of business owners and building permits.
As we mentioned earlier sparsity of data specifically in areas related to spatial analysis makes the process of finding a suitable model challenging. We have investigated and are investigating different spatial matrices to improve the Moran's I criterion. We try to compensate the lack of availability of data set by considering suitable and more statistically informative priors to predict θ ts il and its corresponding Bayesian risk via Gibbs samplers. The best predictors can be obtained from deviance information criteria (DIC) and conditional predictive ordinate (CPO). We provide residual analysis to ensure we have selected the best model and the results would be ready for the JSM 2017.