TY - JOUR
AU - Savitsky, Terrance D
AB - Abstract The US Bureau of Labor Statistics use monthly, by-state employment totals from the Current Population Survey (CPS) as a key input to develop employment estimates for counties within the states. The volatility of monthly CPS state totals, however, compromises the accuracy of resulting county estimates. Typically employed models for small area estimation produce denoised, state-level employment estimates by borrowing information over the survey months, but assume independence among the collection of by-state time series, which is typically violated due to similarities in their underlying economies. We construct Gaussian process and Gaussian Markov random field alternative functional prior specifications, each in a mixture of multivariate Gaussian distributions with a Dirichlet process (DP) mixing measure over the parameters of their covariance or precision matrices. Our DP mixture of functions models allow the data to simultaneously estimate a dependence among the months and between states. A feature of our models is that those functions assigned to the same cluster are drawn from a distribution with the same covariance parameters, so that they are similar, but do not have to be identical. We compare the performances of our two alternatives on synthetic data and apply them to recover denoised, by-state CPS employment totals for data from 2000–2013. 1. INTRODUCTION The Current Population Survey (CPS) is a household survey administered by the Census Bureau for the Bureau of Labor Statistics (BLS) that collects information on the employment status of the civilian noninstitutional population over sixteen years of age residing in the fifty states plus the District of Columbia. The CPS includes sixteen questions, which are administered in-person by a trained interviewer, to ascertain the nature of employment for each member of the household, under a panel survey design where an adult member of each household is interviewed in four successive months, then not interviewed for the following eight months, and included, again, for the following four months. In April 2017, there were 74,248 households in the sample 61,650 of these were eligible households, and of those, 52,996 were actually interviewed. This represents a response rate of approximately 86 percent. The CPS publishes monthly, state-indexed survey direct total employment estimates. The state total employment survey direct estimates are constructed from the sampled household responses in that state by using sampling weights that account for the complex sampling design under which the households were sampled. Our data are these time-indexed survey direct total employment estimates for each state. These direct estimates are volatile, so BLS applies independently specified regression models for each state to reduce the volatility before apportioning these estimates to labor market area or county as part of their local area unemployment survey (LAUS) program. This approach does not exploit the dependencies among the time series. State administrators frequently express concern that the resulting county-level estimates composed from the state-level regression models are overly volatile, in that they believe much of the month-over-month movement is noise, and not related to underlying economic drivers. The LAUS program seeks a more accurate (and less volatile) denoised time series for each state, by exploiting the dependence of CPS direct estimates across states. Bayesian nonparametric Gaussian processes (Vanhatalo, Riihimaki, Hartikainen, Jylanki, Tolvanen, et al. 2014) treat the collection of time series as sums of latent functions and a noise process. The model is focused on latent functions generated from a multivariate Gaussian distribution under covariance constructions that regulate smoothness properties of the resulting functions (Rasmusen and Williams 2006). Vanhatalo et al. (2014) specify a Gaussian process prior that is independent over the states. Dependence among a set of noisy functions over the states can be modeled by combining the functions into a vector of length equal to the number of states multiplied by the number of time-indexed measurement waves. Sampling the full joint posterior distribution would be computationally prohibitive, as computation time scales with order equal of the cube of the resulting vector length. Modeling the functions as a single vector is also not natural, because the ordering of the states may notably impact both global and local smoothness properties of the resulting vector, and this approach assumes a continuous transition between state functions. We offer an approach that builds on the independent nonparametric Bayesian functional models by probabilistically tying the set of latent functions together under a marginal nonparametric mixture prior, where the mixing measure is assigned a Dirichlet process (DP) prior (Ferguson 1973). Our modeling framework is similar to Gelfand, Kottas, and MacEachern (2005), but their approach imposes the Dirichlet process prior directly on the functions, so that co-clustered functions must be exactly the same. This is too restrictive for our application. Instead, impose the DP prior on the covariance matrix parameters of the generating Gaussian distribution, such that state functions which are co-clustered are drawn from a Gaussian distribution with the same parameters, but are not required to be exactly equal. Our proposed approach allows information in the data to extract signal from noise in the time series for each state. We formulate our nonparametric mixture approach using two alternative prior constructions to parametrize the covariance matrix of a multivariate Gaussian prior imposed on the set of latent functions; first, a Gaussian process prior; and second, an intrinsic Gaussian Markov random field prior. These alternative formulations are developed and their properties contrasted in section 2. Schemes to sample the posterior distributions under each of these alternatives are sketched and discussed in section 3. We compare their performance, both in fitting the functions and uncovering a clustering structure used to generate the functions, in section 4. Our nonparametric mixture models are applied for estimation of state-level employment totals obtained from the CPS in section 5, followed by a concluding discussion in section 6. 2. MODELS FOR FUNCTIONAL DATA 2.1 CPS Employment Count Data Our data is composed of T=158 months of (sampling weighted) direct estimates of employment totals for each of N = 51 states (including the District of Columbia) based on the Current Population Survey. The employment totals are influenced by the underlying economic conditions of the states. Performing inference on clusters of distributions generating the set of state-indexed functions provides a context for reported employment and unemployment levels, as we show later. We denote the set of (T=158)×1 vectors of survey direct estimates for a collection of N = 51 states as {yi}i=1,…,N ⁠. We estimate the dependence through a prior construction that permits a grouping or clustering of covariance (under a Gaussian process [GP] prior specification) or precision (under an intrinsic Gaussian Markov random field [iGMRF] prior specification) parameters, that we index by state. These state-indexed covariance or precision parameters are used to estimate denoised, smooth functions {fi} ⁠, from the {yi} ⁠. We first outline GP and iGMRF formulations with simpler, global covariance and precision parameters, respectively, (θ,κ) ⁠, not indexed by state, to illustrate the properties of functions generated under these prior specifications. Employing a single set of covariance or precision parameters assumes the functions are drawn exchangeably. We subsequently extend both of these models using nonparametric mixture formulations over the states for {θi} and {κi} under the GP and iGMRF models, respectively, to allow for dependence among the functions. 2.2 Gaussian Process Model Our Gaussian process (GP) model is parameterized through the covariance matrix, C(θ) ⁠, where θ are estimated by the data and control the trend, length scale (or wavelength), and variation in generated functions, {fi} ⁠. We specify a GP formulation in the following probability model: yiT×1∼indNT(fi,τɛ−1IT), i=1,…,N(1a) fi∼iidNT(0,C(θ)T×T)(1b) θ1,…,θP|a,b∼iidGa(a=1.0,b=1.0)(1c) τɛ|c,d∼iidGa(c=1.0,d=1.0),(1d) where θ=(θ1,…,θP) ⁠. The GP model specifies a formula for the covariance matrix, C(θ) ⁠, with C(θ)=(Cfj,fℓ)j,ℓ∈(1,…,T) Cfj,fℓ=1θ1(1+(tj−tℓ)2θ2θ3)−θ3, for P = 3. This is known as the rational quadratic formula, which may be derived as a mixture over the length scale parameter of more commonly used squared exponential kernels, 1/θ1 exp ⁡((tj−tℓ)2/θ) ⁠, with the inverse of the length scale parameter, θ−1 ⁠, which controls function periodicity, distributed as a Beta distribution with hyperparameters (θ3,θ2−1) (Rasmusen and Williams 2006). The vertical magnitude of a function drawn under the rational quadratic covariance formula is controlled by θ1, while θ2 controls the mean length scale or period of the function, and θ3 controls smooth deviations from θ2. As θ3↑∞ ⁠, this formulation converges to a single squared exponential kernel with length-scale, θ2. So the T × T covariance matrix C(θ) is parameterized through the rational quadratic covariance formula, specified with parameters, θ=(θ1,θ2,θ3) ⁠. Figure 1 displays a randomly generated function under each of 3 distinct value settings for θ=(θ1,θ2,θ3) to provide insight into how θ produces distinct features in f. Our choice of the rational quadratic covariance formulation is intended as a parsimonious specification of a covariance matrix that is able to discover multiple length scales or wavelengths in a function, in lieu of using a set of squared exponential covariance terms (e.g., in an additive formulation), with each term specialized to a particular length scale. For our CPS data, θ2 and θ3 captures economic shocks that are important to capture as features rather than noise. Figure 1. Open in new tabDownload slide RandomlyDrawn Functions, f, from a GP Using a Rational Quadratic Covariance Function Under Distinct Values for θ=(θ1,θ2,θ3) ⁠. The top panel employs θ=(0.4,0.2,5.0) to generate a “rough” surface, while the bottom panel specifies (0.4,4.0,5.0) to characterize a “smooth” surface, and the middle panel employs (0.4,0.2,0.3) to produce a rough surface characterized by a high-degree of length-scale variation. The rational quadratic formula also produces smooth surfaces, {fi} ⁠, that are constrained to be differentiable at all orders and separate signal from noise. The GP of (1) under the rational quadratic covariance formula is very flexible, allowing the estimation of functions, {fi} ⁠, to be non-linear of any order within the space of smooth functions. This GP is equivalent to an infinite dimension polynomial spline, even though the parameterization of the GP is through the covariance matrix, rather than the mean. A linear mean model may be re-parameterized as a zero mean GP by marginalizing over the regression coefficients, which demonstrates that linear surfaces are a subset of the space of smooth functions parameterized in (1) (Neal 2000b). The {yi} are assumed to conditionally independent, given {fi} ⁠. Standard errors for the monthly CPS direct estimates were not available to us. We treat the precision parameter, τɛ ⁠, as unknown and assign it a Gamma prior in our model, which remains fully identified. 2.3 Intrinsic Gaussian Markov Random Field Model An iGMRF prior is typically assigned to the first difference, Δfij=fi,j+1−fi,j∼iidN(0,κ−1) and may be viewed as a probabilistic local smoother indexed by time. Here κ is a precision parameter, analogous to the θ1 precision parameter under the GP formulation. This approach uses nearest neighbors, (fi,j−1,fi,j+1) ⁠, to encode the time-indexed dependence structure among the fij,j=1,…,T for each domain, i∈(1,…,N) ⁠. Modeling first differences may produce an excessively rough (though continuous) surface that over-fits the data, so we employ a prior construction based on the second difference approximation to the first derivative, Δ2fij=Δ(Δfij)∼iidN(0,κ−1) ⁠. This prior on second differences produces the joint distribution, fi∼iidκT−22 exp ⁡(−κ2∑j=1T−2(fij−2fi,j+1+fi,j+2))=κT−22 exp ⁡(−12fi′Rfi), where R=κQ specifies a band-diagonal precision matrix with non-zero entries, (Qj,j−2,Qj,j−1,Qj,j,Qj,j+1,Qj,j+2)=(1,−4,6,−4,1) for 2<j<T−2 for non-boundary parameters (Rue and Held 2005). Under this construction, R is rank-deficient (of rank T – 2) as the rows sum to 0, so that the joint distribution is improper; in particular, the prior for the T×1, fi ⁠, is invariant to the addition of any second order polynomial because the prior supplies no information about such polynomials. We may view the joint distribution as the product of a proper distribution on the space of T−2 differences (by employing the Moore-Penrose pseudo inverse, R− ⁠, and |R| as the product of the T – 2 non-zero eigenvalues of – Q) and an improper, noninformative prior on the order 2 polynomials. These Gaussian Markov random field priors specified through a precision matrix have the property that fij⫫fik|fi,−jk↔Rjk=0 ⁠, which allows for a parsimonious Q, from which we specify a proper set of full conditional distributions that we use for posterior sampling, fij|fi,−j,κ∼N(−1Qjj∑k:k∼jQjkfik,(κQjj)−1)(2a) =N(46(fi,j+1+fi,j−1)−16(fi,j+2+fi,j−2),(6κ)−1),(2b) where {k:k∼j} denotes the set of neighboring time points, {k}, of time, j. The prior mean for each fij is composed as a weighted average of its order 2 nearest neighbors. Rue and Held (2005) refer to this construction as a random walk prior of order 2 or RW2(κ) ⁠. One may equivalently parameterize, R=κ(D−ρΩ) ⁠, for diagonal weight matrix, D, with (diagonal) entries equal to those of Q. Element, j, of D, counts the number of neighbors of time, j, where the function value at a time point with more neighbors will have a relatively higher precision. The adjacency matrix, Ω ⁠, is specified with zeroes on the diagonal, and off-diagonal elements are filled with the negative of values from the off-diagonal elements of Q and encode adjacency relationships among the time points (Banerjee, Wall, and Carlin 2003). Under this parameterization, ρ, −1<ρ<1 ⁠, is the autocorrelation parameter that shrinks the mean to 0, but produces a proper joint prior for fi ⁠. The RW2(κ) prior construction, however, sets ρ = 1. As with (1), we specify τɛ∼Ga(c=1.0,d=1.0) ⁠. The GP model (1) includes parameters (θ2,θ3) for the length-scale, implying more flexibility than the iGMRF model, which hard codes the length scale in the prior for Q. We assess this question in section 4. 2.4 Accounting for Dependence among Functions We introduce an extension of (1), which indexes the GP covariance formula parameters, {θi}={(θi1,…,θiP)} ⁠, by state, i∈(1,…,N) ⁠, to permit their probabilistic clustering. Under our rational quadratic formulation with P = 3, each θi=(θi1,θi2,θi3) ⁠, so that each state receives its own GP covariance formula parameters (which parameterizes the T × T covariance matrix of the GP prior for domain, i). These sets of parameters are clustered over the states under the following specification, fi∼NT(0,C(θi))(3a) θ1,…,θN|G∼iidG(3b) G∼DP(α,G0),(3c) where {θi}i=1,…,N receive a random distribution prior, G, drawn from a Dirichlet process (DP), with a concentration parameter, α, that controls the amount of variation in G around prior mean, G0. The base or mean distribution, G0, is typically constructed as parametric; in our case, G0=∏p=1PGa(1.0,1.0) ⁠, which we used under the global GP model of (1) parameterizing a single vector, θ ⁠, for all states. Equation 3 describes a mixture model of the form, f|G∼iid∫NT(0,C(θ))G(dθ) ⁠, where G is the mixing measure over the P×1 covariance parameters, θ ⁠. We examine the clustering property of the DP by expressing it in the discrete, stick-breaking form of Sethuraman (1994), G=∑h=1∞phδθh∗,(4) a countably infinite mixture of weighted point masses, where “locations,” θ1∗,…,θM∗ ⁠, index the unique values for the {θi} ⁠, with M≤N (states) that we interpret as clusters. The maximum number of clusters assigns each state to its own cluster, which countably increases with N. We record state cluster memberships with s=(s1,…,sN) where si=ℓ denotes θi=θℓ∗ so that (s,{θm∗}) provides an equivalent parameterization to {θi} and we recover θi=θsi∗ ⁠. We conduct posterior sampling with the cluster locations and assignments, rather than directly sampling {θi} ⁠, as the former produces notably better mixing because it separates re-assignments to clusters from updates to the values. The weight, ph∈(0,1) ⁠, is composed as ph=vh∏k=1h−1(1−vk) where vh is drawn from the beta distribution, Be(1,α) ⁠. This construction provides a prior penalty on the number of mixture components. A higher value for α, however, will generate smaller values for {vh} ⁠, and hence, “breaks” off more clusters (with unique locations), {ph} ⁠, from the unit stick. We place a further prior, α∼G(1,1) ⁠, to allow posterior updating in recognition of the relatively strong influence it conveys on the number of clusters formed (Escobar and West 1995). The DP construction assumes exchangeability of the (θi) ⁠, a priori, given random measure, G, but the almost surely discrete construction of the DP produces estimates that are not exchangeable, a posteriori. The prior specification for cluster assignments, (si|s−i) ⁠, (under our re-parameterization to {s,(θm∗)m=1,…,M} achieved by marginalizing over G), induces a uniform probability for co-clustering among the states. Let Cτ(θm∗)=C(θm∗)+(1/τɛ)IT after integrating out fi from (1a) to produce the marginal likelihood, yi|s,(θm∗)∼NT(0,Cτ(θsi∗)) ⁠. If the likelihood values for two states, j and k, NT(yj|0,Cτ(θs∗)) and NT(yk|0,Cτ(θs∗)) ⁠, are both relatively high for assignment to cluster, sj=sk=s ⁠, due to underlying similarities in their economies or due to other factors related to their closeness in their geographic locations, then the posterior probability for co-clustering states j and k will be relatively high (versus uniform, a priori). This posterior estimation mechanism conducts unsupervised (probabilistic) clustering. Sharper (lower posterior variance) estimates over the space of clusterings may be obtained, particularly when the number of observations are higher (than our N = 51), by indexing either the weights or locations in (4) to include predictors in lieu of our unsupervised formulation. An analogous extension of (2) is fij|fi,−j,{κi}i=1,…,N∼N(−1Qjj∑k:k∼jQjkfik,(κiQjj)−1)(5a) κ1,…,κN|G∼iidG(5b) G∼DP(α,G0),(5c) which may be expressed as fij|fi,−j,G∼iid∫N(−1Qjj∑k:k∼jQjkfik,(κQjj)−1)G(dκ) ⁠, marginalizing over {κi}i=1,…,N ⁠. We specify the base distribution G0=Ga(1,0.1) ⁠, (which generates locations, {κm∗} ⁠) as weakly informative with a large variance. We select these hyperparameter settings because they produce a larger prior variance for κ in (5) than for Θ under (3). The parameterization of the precision matrix of each function using only a single parameter (κ) makes the posterior sensitive to the prior specification. As with (3), we sample κ indirectly through cluster assignments, s, for the N×1 states, and location values κ1∗,…,κM∗ ⁠. We also specify the same prior for the DP concentration parameter, α∼Ga(1.0,1.0) ⁠. The stick-breaking formulation illustrated in (4) to induce the unknown measure, G, may be easily generalized beyond the DP by varying the form for the Beta prior assigned to vh; for example, the Poisson-Dirichlet (PD) alternative of Pitman and Yor (1997) includes an additional “discount” parameter, δ, along with the DP concentration parameter, α, to specify the hyperparameters of the Beta distribution for vh that tends to generate more locations, a priori, inducing a less informative prior for the number of clusters (which is random under DP and PD constructions). Large clusters, however, receive more prior weight under the PD than for the DP. The DP is a special case of the PD, and both may be located within a larger class of generalized gamma processes (Lijoi, Mena, and Prünster 2007). The PD and DP produce nearly identical posterior results under our data application due to a high level of information in the data, so we focused our exposition on the simpler model. 3. COMPUTATION The mixtures of Gaussian processes formulation of (3) is far more computationally intensive than the mixtures of iGMRFs model presented in (5) because computing the cholesky decomposition of the T × T covariance matrix is O(T3) ⁠. We report computation time comparisons in section 4, where we also demonstrate that the mixtures of GPs is more successful at discovering the correct number of generating clusters than is the mixtures of iGMRFs. So we are motivated to mitigate the computational burden associated to drawing posterior samples under a GP prior on the functions. A typically used approach to improve computation for a GP probability model employs a sparse approximation to the GP that utilizes some subset of the time points as “inducing” inputs (Vanhatalo et al. 2014). We expect functions, {fi} ⁠, for our CPS data to each express multiple length scales because our estimation period encompasses over one hundred months, such that a sparse approximation may fail to capture key local features, such as the economic shocks we earlier mentioned. So we adapt a recently developed posterior sampling algorithm of Wang and Neal (2013) to our mixtures of GPs, which generates a sequence of proposals for locations, {θm∗} ⁠, in a lower dimensional space using a subset of the T time points to build the GP covariance matrix, C. This reduces computation time (by five- to 10-fold versus sampling the full-dimensional covariance matrix). The proposed moves in the lower dimensional space may be sampled in a rapid fashion before stepping back up to the full space to compose a probability-of-move under a Metropolis-Hastings sampling scheme. Yet, sample draws under this method will be from the exact posterior distribution, rather than from a sparse approximation, to capture local-in-time features. Appendix A provides an exposition of the posterior sampling algorithms used for both the mixtures of GPs and iGMRFs models, which we have implemented in the growfunctions package for R Core Team (2014), using C++ to speed computation. The package is fully documented and available on CRAN for downloading and includes the data used in this paper. 4. SIMULATION STUDY We conduct a simulation study with the goals to: 1). Assess and compare the accuracy of estimating denoised functions between our GP and iGMRF formulations; 2). Assess and compare the accuracy of detecting clusters or sub-groups of denoised functions. We generate each latent synthetic function, f (from which we render the observed time series, y) under different parameterizations than either of our models as a means of assessing the relative adaptability and robustness of our estimation models to real-world data generating processes. Our first simulation focuses on generating relatively complicated surfaces from a mixture of K = 2 length scales specific to each of M = 3 clusters. We randomly (and disjointly) allocate N = 100 domains to M = 3 clusters for T = 158 time points and generate denoised functions, {fi} ⁠, in each cluster using the addition of two squared-exponential terms, fi∼NT(0,C(θ1,si∗)+C(θ2,si∗)), where the (⁠ j,ℓ ⁠)th cell in the T × T matrix C is Cfij,fiℓ(θk,si∗)=1θk,si,1∗exp ⁡((tij−tiℓ)2θk,si,2∗), The indicator k∈{1,2} denotes the covariance term with associated locations, {θk,m,1∗,θk,m,2∗} ⁠, for a squared exponential covariance (which has two location parameters). Label (k,m,1) indicates a location parameter that controls the vertical scale for covariance term k in cluster m, and (k,m,2) denotes a location parameter that controls the length scale or degree of smoothness. The vector s=(s1,…,sN), si∈(1,…,M) ⁠, assigns domains to clusters, so that the parameter controlling the vertical scale in covariance term, k = 1, for domain i is θ1,si,1∗ ⁠. The columns of the matrix, θ∗       m=1  2        3Θ∗=θ1,m,1∗θ1,m,2∗θ2,m,1∗θ2,m,2∗[2.610.380.913.003.531.561.042.260.840.220.150.71] specify the set of four locations generated for each of the M = 3 clusters that we used to construct {fi} ⁠. Rows θ1,m,2∗ and θ2,m,2∗ of this matrix imply that the surfaces in each cluster are formulated from the sum of a long length scale, or smooth term, and a short length scale, or rough term. So each surface is drawn from a mixture of two scales. The long length scale locations were generated from, θ·,·,2∗∼Ga(3,2) ⁠, while the short length scale locations were drawn from, θ·,·,2∗∼Ga(2,5) ⁠. The vertical scale parameters for both terms were drawn from, θ·,·,1∗∼Ga(3,3) ⁠. The observed surface for domain i is modeled as yi∼NT(fi,τɛ−1IT) ⁠, where the global noise precision parameter, τɛ ⁠, is set to produce a 20 percent noise-to-signal ratio (based on the average variance of the {fi} ⁠). The first of our two primary inferential goals is to uncover the denoised latent functions, {fi} ⁠. We randomly set 10 percent of the NT observations to missing to use them as a test set. The GP and iGMRF models provide estimates, {f^ij} ⁠, for each of these missing values, and we compute a mean-squared prediction error (MSPE) based on the squared difference between the estimated and true function values for the test set. We then normalize the MSPE by the variance of the test set latent functions to produce a normalized MSPE ∈(0,1) ⁠, where larger values indicate a relatively higher error rate. Figure 2 compares fitted results (solid lines) to data values (circles) for the iGMRF RW2(κ) (“random walk of order 2”) model, in the left column, with the GP formulation under a single rational quadratic covariance formula, in the right column. The rows represent each of the M = 3 clusters, and the plot in each cell is for a randomly selected domain from the represented model-cluster combination. We see that the GP model tends to attenuate the peaks to a greater degree than does the iGMRF. The normalized MSPE values are 0.39 for the GP and 0.54 for the iGMRF. The GP fits notably better, because it avoids over-fitting. The GP model employs three parameters that regulate the estimated functions, while the iGMRF only employs a single parameter; in particular, the GP covariance formula allows the data to estimate the length scale and local deviations from that scale, while the length scale is hard-coded or set in the (order 2 or RW2) iGMRF model. Figure 2. Open in new tabDownload slide Posterior Mean Values for Estimated Functions, {fi}, from the iGMRF Model (“rw2” for RW2), in the Left-Hand Column, and Under the GP (“rq” for Rational Quadratic Covariance Formula), in the Right-Hand Column, under a Simulated Data Set Generated from a GP Mixture of Two Squared Exponential Covariance Terms. The rows represent each of the M = 3 clusters and each plot includes a randomly selected domain. Our second inferential goal is to utilize the clustering to provide context about how the functions compare among the observed domains (e.g., geographic or industry labels). Both the GP and iGMRF models employ a DP mixture formulation under which each posterior sampling iteration assigns the covariance or precision parameters to clusters. The functions, fi ⁠, for those domains assigned to the same cluster are drawn from a GP with the same covariance formula parameters. So the domain-indexed functions are not identical, but are similar, because they are drawn from the same multivariate Gaussian (GP) distribution. Both the number of clusters, M, and the allocations of domains to those clusters may change with each posterior sampling iteration, so that these samples, taken together, formulate draws from the posterior distribution over the space of clusterings (or partitions). The relative concentration of this distribution may be assessed by examining the degree of similarity in number of clusters and assignment of the domains to those clusters among the posterior draws. We select one clustering of domains from among the MCMC draws by using the least-squares algorithm of Dahl (2006), which builds an N × N matrix of pairwise clustering probabilities to summarize the posterior draws and selects the clustering that is “closest” to this matrix under a squared Euclidean distance metric. Once we have selected a clustering, we build an N × N pairwise matrix, where cell (i, j) is filled with a 1 if states i and j are clustered together, and 0 otherwise. We also build this pairwise clustering matrix for the true clustering used to generate the observed (synthetic) values. We then compare the true and model-based pairwise matrices and compute a percent mis-clustering based on the sum of the cell values that don’t agree. Figure 3 presents a visual schematic of the true clustering for our synthetic data and the estimated clustering structure under each of the two models. Each plot panel (in each of three groupings of panels) collects the posterior mean values for estimated functions, (fi) ⁠, assigned to the cluster for that panel. We recall that co-clustered domains have similar but not identical estimated functions. Each solid line in a panel represents the posterior mean value for a function (over the T = 158 time points) associated to a domain. The collection of solid lines in a panel display the (posterior mean values for) functions of the collection of domains assigned to the cluster indexed by that panel. The dashed line in the panel is the mean of the collection of posterior mean function values for those domains assigned to the cluster. So the dashed line provides an indication of a common trend among the co-clustered functions (and it is not a posterior mean, but the mean over a domain-indexed collection of posterior mean function values). Each grouping of panels, together, represents a clustering (composed as a set of clusters) under a different model. The top row presents the true clustering under three clusters used to generate the observed {yi} ⁠, while the second row presents the estimated clustering for the GP model (selected under the least squares algorithm), and the last two rows present the clustering for the iGMRF model. We see that the GP model reproduces the generating cluster well, with a mis-clustering error of only 2 percent. The iGMRF model discovers the three clusters, but also creates echoes of the third cluster, producing a total of five clusters and a resulting mis-clustering error of 20 percent. Perhaps even more important than the correct estimation for the number of true clusters, M, is the accuracy with which pairs of domains are correctly co-clustered; because if domains are relatively correctly clustered with one another, then the inference about clusters is more likely to be useful, even if a true cluster is split between two clusters in the estimation model. Our mis-clustering rate, however, indicates that the iGMRF, while still correctly clustering together 80 percent of domains, yields substantially more error than does the GP. Figure 3. Open in new tabDownload slide Cluster Sub-Groupings of the Posterior Means for Estimated Functions, {fi} ⁠. Each panel plots the subset of functions assigned to a cluster under a simulated data set generated from a GP mixture of two squared exponential covariance terms. Each line in the panel is a posterior mean value for a domain-indexed function and all functions associated with domains assigned to the cluster indexed by the panel are co-plotted. A loess smoother is added in each panel with a dotted line to highlight patterns among the assigned functions. Each row represents the clustering or partition under a model. The first row presents the true functions and clustering, while the second row presents the estimated clustering from GP model and the last two rows, the GMRF model. We expect the GP model with a rational quadratic covariance formulation to perform relatively better than the iGMRF model since the rational quadratic is designed for estimation under varying scales. Our next simulated procedure generates latent functions from a model closer to the iGMRF. Recall that the T × T precision matrix of GMRF may be expressed as, Ri=κsi∗(D−ρΩ) ⁠, where we set the strength of correlation parameter, ρ = 1 (which produces a rank-deficient precision matrix), to produce our iGMRF formulation. If we set 0<ρ<1 ⁠, we recover a full rank precision matrix, but at the expense of shrinking the conditional mean values to zero. So we next draw our true function from, fi∼NT(0,κsi∗(D−0.95Ω)) to construct our next synthetic dataset, under a second order specification for {D,Ω} ⁠, as described in section 2.3 using the same N, T, M as for the first simulation. Locations, {κm∗}m=1,…,3 ⁠, are generated from a Ga(1,1) ⁠. We again set τɛ to produce a 20 percent noise-to-signal ratio. The resulting functions are similar to those that would be generated under an iGMRF. While the iGMRF achieves an nMSPE of 0.159, much lower than in the first simulation scenario, the GP achieves a similar level of 0.157. The similar results suggest robustness of the GP model to alternative true function values. We next estimate the GP formulation of (1), where we exclude the modeling of clusters among {fi} to assess its relevance. We compare MSPE fit performances between including and excluding the clustering prior on the second simulated dataset. The resulting normalized MSPE of 0.3 when excluding the clustering prior is notably higher than the value for both the DP mixtures of GPs (with an MSPE of 0.157) and GMRFs (with an MSPE of 0.159) of (3) and (5). So removing the ability for the data to estimate sub-groups among the generating distributions for the domain-indexed functions (e.g., states) reduces estimation accuracy because all of the states are assumed to express the same underlying patterns in their true functions. We conclude our simulation study with an exploration of the repeated sampling (frequentist) properties of our two Bayesian estimators for the mixtures of GPs and iGMRFs formulations. Function values, (fi)i=1,…,(N=300) ⁠, are generated using a single covariance term under a squared exponential covariance formula with each domain, i∈(1,…,(N=300)) ⁠, randomly allocated to 1 of M = 3 clusters with associated cluster locations given by,                 θ∗ m=1  2    3Θ∗=θm,1∗θm,2∗[1.502.001.03.004.001.0]. The right-most cluster (covariance formula) location parameters will induce relatively rough or high-frequency functions, while the location parameters for the next cluster from the right will generate relatively smooth functions and the left-most cluster location parameters will generate generate surfaces of smoothness in between the first two. Each fi=(fi1,…,fiT) is synthesized for T = 40 months. We next sample the associated yi|fi,τɛ as before, setting the global precision parameter, τɛ ⁠, to produce a 10 percent noise-to-signal ratio. We conduct B = 100 Monte Carlo simulation runs, where a new data set, yb,i, b=1,…,B, i=1,…,N ⁠, is sampled on each Monte Carlo iteration from the same fi ⁠, such that the same true function values are employed across all of the Monte Carlo iterations. We assess the relative fit performances of the iGMRF and GP models by randomly excluding 10 percent of the NT observations. The same pattern of missingness is employed across the Monte Carlo simulations. We assess the performance of the models by seeing how well they predict the true function values, (fij) ⁠, for those months, j, where the associated data values are treated as missing (and excluded from model estimations). Each panel in figure 4 displays the mean across the B = 100 Monte Carlo iterations of posterior mean estimated function values under the GP (gp_se) and iGMRF (gmrf_rw2) as compared to the true function values (true), fi ⁠, for randomly selected domains, i, in a cluster, m∈(1,…,(M=3)) ⁠. We include the full uncertainty of posterior estimation and repeated sampling by incorporating the 95 percent credible intervals formed by concatenating the posterior draws for the estimated functions over all of the B = 100 Monte Carlo iterations. Figure 4. Open in new tabDownload slide Credible Intervals of the Posterior Draws for fiOver the B = 100 Monte Carlo Iterations for a Single Domain in Each of M = 3 Clusters. Each plot panel displays the 95 percent credible intervals constructed by concatenating the posterior draws over the B = 100 Monte Carlo iterations for a randomly selected fi in each of M = 3 clusters. Each panel displays the true function value, fi ⁠, labeled as “true,” the estimated function value and associated credible intervals from the GP model (gp_se) and the same for the iGMRF (gmrf_rw2) model. This figure reveals that, under repeated sampling, both models tend to over-smooth in predicting the true function values for the excluded months in each domain. The GP model tends to predict the true values better, displaying less smoothness for the predicted (missing) values. We see that the GP intervals are generally shorter and cover the truth better than the iGMRF intervals, though both do well in smoother portions of the true functions. The credible intervals for the GP and iGMRF notably depart in portions of each true function which are less smooth/rougher, and the GP better covers the true values in the case of greater roughness in the true function. This is expected, because the GP employs a length-scale parameter, which is not part of the iGMRF model. Yet, as seen in the previous simulations, GP tends to produce smoother estimates of fij for the observed months, j, than iGMRF. The higher degree of smoothing of GP on observed values induces more flexibility in the fit on the predicted values by adaptively discovering longer or shorter memory patterns in the underlying signal. Figures 5 and 6 provide information on the relative variation of the fit performances over the B = 100 Monte Carlo simulations, by displaying the distributions over repeated sampling of the normalized MSPE (for the predicted months) and the misclassification rate, respectively. The superior fit of GP is partly achieved through better clustering accuracy. Figure 5. Open in new tabDownload slide Repeated Sampling Distribution of the Normalized MSPE Statistic Over the B = 100 Monte Carlo Simulations. Figure 6. Open in new tabDownload slide Repeated Sampling Distribution of Clustering Misclassification Rate Statistic Over the B = 100 Monte Carlo Simulations. The GP model has higher computational cost than the iGMRF model. We roughly achieve the same information (effective sample size) in 10,000 iterations of the GP model as we do in 25,000 iterations of the iGMRF; the former consumes 28,900 CPU-seconds and the latter 1,226 CPU-seconds, on a single core of an Intel-i7–3450-powered laptop, with the estimation routines both written in C++. As mentioned, our tempered posterior sampling approach for the GP makes proposals for the covariance formula parameters in a lower dimensional space; this reduces computation time by an order of magnitude of five. The successively coarser covariance approximations for likelihood evaluations in the temporary space were computed with (100, 60) of the T = 158 time points. It is likely we could do better if we were to optimize these choices. While sparse approximations are used to improve the computational scalability of the GP, local-in-time features, such as economic shocks in our CPS appplication dataset, may be lost or overly smoothed. One has to evaluate the utility of a sparse method for each dataset, while our faster sampling scheme may be employed with any dataset. The computation time for the GP, while relatively high, does fit within BLS monthly production schedules. 5. CPS EMPLOYMENT APPLICATION We return to our motivating application, the T = 158 month series of employment totals for N = 51 states (including the District of Columbia) from the CPS. Each state’s observed employment series is standardized in order to facilitate comparisons across states based on similarities in the shapes and patterns expressed. We perform inference on the distribution over clusterings to understand which states express similarities in their employment patterns. We fit the GP mixture model of (3) to the CPS data. Figure 7 presents the allocation of posterior mean values for the estimated latent functions, {fi} ⁠, for states, i=1,…,(N=51) ⁠, into assigned clusters of the clustering selected under the least squares criterion of Dahl (2006). We recall that the least squares algorithm selects a single clustering from among the set of posterior samples for the clustering. Each panel of the figure displays the posterior mean values for the estimated, denoised functions for those states assigned to each of two clusters estimated by the model. The dashed line in each panel averages over the collection of state-indexed posterior mean function values in each month and represents an average trend expressed by the functions of co-clustered states. The dashed lines in each panel express less variation than the collection of state-indexed posterior mean function values due to the averaging. Examining the pattern in the collection of functions assigned to each cluster reveals that the clusters are primarily differentiated by the employment sensitivities expressed by member states to the “great recession” (a local-in-time phenomenon), that began in approximately 2008 (at month 96). The latent functions for those states included in the left-hand panel express a notably longer-lasting trough than those states included in the right-hand panel (where panels index cluster memberships). The initial rate of employment drop appears similar in both clusters, but the decline continues longer and recovers more slowly among states assigned to the left-hand cluster. We conducted additional runs that varied the shape hyperparameter of the DP concentration prior on α∼Ga(a4,1) ⁠, with a4∈{1−4} to induce successively higher prior number of clusters, but the estimated GP models over the differing values of a4 all selected the same M = 2 number of clusters and assignment of states to those clusters. Estimation under the iGMRF mixture model of (5) produces the same selected clustering, including the number of clusters and assignment of states to clusters, which indicates the data strongly informs the estimated joint posterior distribution over the space of partitions. Figure 7. Open in new tabDownload slide Posterior Mean Estimated Functions, {f}, Grouped into Clusters. Each panel plots the set of functions assigned to a cluster for the CPS state-level dataset under the selected least squares clustering. A loess smoother is added in each panel with a dotted line to highlight a pattern among the assigned functions. The left-hand cluster includes a number of states that experienced so-called bubbles in property values, such as California, Florida, Hawaii, and Nevada, as well as states sensitive to both financial services and auto industries, such as New York and Michigan. (See US Bureau of Labor Statistics [2012] for more background and discussion). Table 1 lists the thirty-eight states assigned to the left-hand cluster of figure 7 and the thirteen states assigned to the right-hand cluster. Each panel in figure 8 presents the fitted latent functions from the GP model of (3) with a solid line and the associated credible intervals shaded in gray, compared to the actual observed data values. The states represented in the top row of panels are assigned to the left-hand cluster of figure 7, while the states presented in the bottom row are assigned to the right-hand cluster. Examining the states assigned to the left-hand cluster, both New York and Michigan experienced a relatively more rapid descent in employment, despite the large volatility in New York, while Florida appears to have experienced a large, nearly sudden, drop followed by a stabilization and recovery. Contrastingly, Arizona and Massachusetts appear to have rather rapidly bounced up from the initial drop, while Ohio experienced a less dramatic fall in employment. The cluster memberships provide state administrators context on the common economic drivers for their employment trends through the great recession. Table 1. State Memberships in 2 Estimated Clusters Cluster . Member states . Left AK, AL, AR, CA, CO, CT, DC, FL, GA, HI, IA, IL, IN, KS, KY, LA, MD, MI, MN, MO, MS, ND, NE, NV, NY, OK, OR, PA, RI, TX, UT, VA, VT, WV, WY Right AZ, DE, ID, MA, MT, NC, NJ, OH, SC, SD, TN, WA, WI Cluster . Member states . Left AK, AL, AR, CA, CO, CT, DC, FL, GA, HI, IA, IL, IN, KS, KY, LA, MD, MI, MN, MO, MS, ND, NE, NV, NY, OK, OR, PA, RI, TX, UT, VA, VT, WV, WY Right AZ, DE, ID, MA, MT, NC, NJ, OH, SC, SD, TN, WA, WI Open in new tab Table 1. State Memberships in 2 Estimated Clusters Cluster . Member states . Left AK, AL, AR, CA, CO, CT, DC, FL, GA, HI, IA, IL, IN, KS, KY, LA, MD, MI, MN, MO, MS, ND, NE, NV, NY, OK, OR, PA, RI, TX, UT, VA, VT, WV, WY Right AZ, DE, ID, MA, MT, NC, NJ, OH, SC, SD, TN, WA, WI Cluster . Member states . Left AK, AL, AR, CA, CO, CT, DC, FL, GA, HI, IA, IL, IN, KS, KY, LA, MD, MI, MN, MO, MS, ND, NE, NV, NY, OK, OR, PA, RI, TX, UT, VA, VT, WV, WY Right AZ, DE, ID, MA, MT, NC, NJ, OH, SC, SD, TN, WA, WI Open in new tab Figure 8. Open in new tabDownload slide Posterior Mean Values for Estimated Functions, {f}, from the Mixtures of GPs Model, Represented by the Solid Lines, with 95 Percent Credible Intervals Displayed as Shaded Region Compared to the Noisy Series, {y}, Represented by the Hollow Circles, Connected with Dashed Lines. Each panel represents a series for a single state. The state plots in the top row are assigned to the left-hand cluster of figure 7, while the bottom row plots are from the right-hand cluster. The solid line in figure 9 displays the denoised estimated function for the state of Nebraska compared to the noisy CPS estimates (displayed in the hollow circles connected with a dashed line) to demonstrate the reduction of noise accomplished under our mixture of GPs model. BLS have received very positive subjective feedback from states regarding the quality of the signal extracted (from the estimated set of T×1 denoised functions, {fi} ⁠). Figure 9. Open in new tabDownload slide Posterior Mean Estimated Function for the State of Nebraska from the Mixtures of GPs Model, Represented by the Solid Line, with 95 Percent Credible Intervals Displayed as Shaded Region Compared to the Noisy Series, {y}, Represented by the Hollow Circles, Connected with Dashed Lines. Lastly, our use of Dahl (2006) to select a single MAP partition or clustering from the posterior distribution over the space of partitions does not express the uncertainty in this distribution. It is challenging to express the uncertainty around the MAP partition as each partition is composed of assignments of states to clusters. We employ the 51 × 51 matrix of pairwise co-clustering probabilities, which are used to select a single partition under Dahl (2006), to illustrate the variation around the MAP partition. Figure 10 plots the histogram of values for pairwise co-clustering probabilities, πij,i,j∈1,…,N=51 ⁠. The extent of the distance or spread between concentrations of low and high co-clustering probabilities provides an indication of the degree of concentration around the MAP point clustering estimate. We observe in figure 11 concentrations of co-clustering probabilities bounded away from 0.5 to either relatively low or high values, which indicates the distribution of partitions is relatively concentrated around our MAP estimate. See Appendix A for additional plots, including a 51 × 51 heat map of the co-clustering probabilities among the states and an MCMC trace plot for the draws of M, the total number of clusters. Figure 10. Open in new tabDownload slide Histogram of Empirical Pairwise Assignment Probabilities, πij∈[0,1], for the States, Estimated from the post-in MCMC Draws of Cluster Assignments from the Posterior Distribution Over Clusterings/Partitions. A larger spread between low and high values for pairwise assignment probabilities indicates lower variability in cluster assignments. Figure 11. Open in new tabDownload slide Empirical Co-Clustering (Pairwise Assignment) Probabilities for States, πij∈[0,1], Estimated from the post-in MCMC Draws of Cluster Assignments from the Posterior Distribution Over Clusterings/Partitions. The focus is on the co-clustering probability for state i with state j, and not the assignment to clusters. 6. DISCUSSION Our task was to estimate latent, state-indexed functions that improve the efficiency and interpretability of county-level employment statistics constructed from the state CPS estimates. We developed nonparametric mixture formulations that simultaneously estimate the latent, state-indexed functions and allow the data to discover a general dependence structure among them that borrows estimation strength to improve precision. Our simulation study results demonstrated that failing to account for a dependence structure among states in the estimation model lessens the ability to uncover the true latent functions. Our DP mixture of GPs or iGMRFs, outlined in (3) and (5), employ an unsupervised approach for discovering dependence among the state functions based on similarities in the time-indexed patterns expressed in the data. They perform well on our CPS employment count application, uncovering a differentiation among states based on their employment sensitivities to the Great Recession. The simulation study revealed a greater robustness, both in estimation of accuracy of the latent functions and their clustering properties, of the GP mixture model as compared to the iGMRF mixture, due to the regulated smoothness properties of the rational quadratic covariance formulation and its inclusion of more parameters than iGMRF to reflect the trend, scale and frequency properties of the estimated latent functions. The iGMRF computes much faster, however, so it may still be useful, particularly in the case where the clusters are differentiated based primarily on vertical magnitude of the functions. The computing time of the GP mixture is reduced by our adaptation of the faster-computing MCMC algorithm of Wang and Neal (2013). The authors wish to thank colleagues at the Bureau of Labor Statistics who supported this project and provided the data for this project. We thank the following important contributors: Sean B. Wilson, Senior Economist, who formulated the project and provided feedback from the states on our results; Garrett T. Schmitt, Senior Economist, who helped us think through alternative approaches; Bradley A. Jensen, Senior Economist, who provided us multiple data slices that allowed us to perform our estimation. APPENDIX A. POSTERIOR SAMPLING ALGORITHMS A.1 Gaussian Process Sequential Scan We sample the parameters of the mixtures of GPs model, λ=(Θ∗,s,α,τɛ) ⁠, in a sequential scan from the full conditionals after marginalizing over the latent, denoised functions, {fi} ⁠, which is easy to do and leads to slightly better mixing. We must, however, co-sample the {fi} in the case of intermittent missingness-at-random in {yi} in order to allow the sampling of missing values from their posterior predictive distribution. We first outline our sequential scan in the case where we marginalize out the latent functions and conclude this section by highlighting changes required in the case we co-sample them. The mixtures of GPs model specification are highly non-conjugate. We carefully configure blocks of parameters to permit application of posterior sampling approaches designed to produce robust chain mixing under non-conjugate specifications. 1. Sample cluster locations, {θpm∗}p=1,…,P; m=1,…,M ⁠, where p denotes parameter type (e.g., (θ1∗,θ2∗,θ3∗) ⁠) and m=1,…,M denotes the cluster and model error precision, τɛ ⁠: We sample the posterior distribution for locations in by-cluster groups, {θpm∗}p=1,…,P ⁠, from the subset of observations (states) assigned to that cluster because θpm∗⫫θpm′∗ for m′≠m ⁠, a posteriori, in a Metropolis-Hastings scheme using the following log-posterior kernel, log  π(θpm∗|θ−pm∗,s,τɛ,{yi:si=m}) ∝−12nm log(|Cτ(θm∗)|)−12∑i∈smyi′Cτ(θm∗)yi+(a−1) log ⁡(θpm∗)−bθpm∗,(6) where Cτ(θm∗)=C(θm∗)+(1/τɛ)IT ⁠, sm collects the states assigned to cluster m, nm denotes the number of states assigned to cluster m, and (a, b) are shape and rate hyperparameters of a gamma prior, respectively. This posterior representation is a relatively straightforward Gaussian kernel of a non-conjugate probability model. We adapt a Metropolis-Hastings algorithm of Wang and Neal (2013) designed to speed computation for sampling each θpm∗ (and also, τɛ ⁠) by introducing a lower-dimensional temporary space where the likelihood (e.g., the T × T, Gaussian process covariance matrix, C) is approximated using a subset of the T time-points. Starting with the previously sampled value, x^0 ⁠, (for θpm∗ or τɛ ⁠) from the full-dimensional or exact space, the sampling algorithm builds a sequence of proposed transitions by first “stepping up” in the temporary space using increasingly coarse, “tempered” transitions, (Z^1,Z^2,…,Z^n) ⁠, that generate computationally fast approximations for C. These approximations use fewer observations in each step; for example, we apply n = 2 transitions in the lower-dimensional space and use 100 of the T = 158 time points to formulate Z^1 and then 60 time points for Z^2 ⁠. This sequence of coarser transitions is followed by transition steps that employ progressively finer distributions (in reverse order), (Zˇ2,Zˇ1) that” step down” to guide the chain back toward the full dimensional space until we conclude the sequence by proposing xˇ0 from Zˇ1 to be evaluated in the full-dimensional space. We use univariate slice sampling of Neal (2000a) to accomplish these (reversible) transitions in the lower-dimensional space. The proposal (reversible sequence of steps) is then accepted with (for n = 2 tempered transitions) probability of move formulation, min⁡(π1(x^0)π2(x^1)π1(xˇ1)π1(x^1)π2(xˇ1)π1(xˇ0)π(xˇ0)π(x^0)), where “x” denotes the proposals, and “π” posterior kernel evaluations, for each θpm∗ with subscript 0 pertaining to the exact space and (1, 2) to the successively coarser transition distributions in the temporary space in the sequential order of application. If our lower dimensional approximations are relatively good, this approach will speed chain convergence by producing draws of lower autocorrelation since each proposal includes a sequence of moves generated in the temporary space. Since the moves in the temporary space are executed with fast approximations, Wang and Neal (2013) show that this algorithm has the potential to substantially reduce computation time, as compared to the usual Metropolis-Hastings algorithm, for drawing an equivalent effective sample size. The probability of move formulation evaluates the proposals on the full space of size T, however, so that the resulting sampled draws are from the exact posterior distribution, rather than from a sparse approximation. Our adaptation configures Wang and Neal (2013) to apply to clusters of Gaussian processes (rather than to just a single GP) by exploiting the between cluster posterior independence of the {θpm∗} ⁠. 2. Sample cluster assignments, s=(s1,…,sN) ⁠: We marginalize over G in Equation 3, that results in the Pólya urn representation of Blackwell and MacQueen (1973), to sample s from its full conditionals, π(si=s|s−i,Θ∗,α,τɛ,yi)∝{n−i,sn−1+αNT(yi|0,Cτ(θs∗)) if 1≤s≤M−α/c*n−1+αNT(yi|0,Cτ(θs∗)) if s=M−+h,(7) where n−i,s=∑i′≠i1(s(i′)=s) is the number of states, excluding state i, assigned to cluster s, so that states are assigned to an existing cluster with probability proportional to its “popularity”. The posterior assigns a state (through si) to a new cluster with probability proportional to αd0 under the mixture prior in the case of a conjugate formulation. The conjugate specification requires the likelihood to be integrable in closed form with respect to the base distribution, G0, to compute d0=∫N(y|θ,…)G0(dθ) ⁠, which is not the case under a (nonconjugate) GP construction. So we employ the auxiliary Gibbs sampler formulation of Neal (2000b) and sample c∗∈N locations from base distribution, G0, ahead of any assigned observations (e.g., not yet linked to any state), to define h=M−+c∗ candidate clusters in an augmented space. We then draw si from this augmented space, where any location not assigned states (over a set of draws for s) is dropped. This procedure effectively performs a Monte Carlo integration of the likelihood with respect to the base distribution. See Savitsky and Vannucci (2010) for a detailed example of a DP implementation on a GP. The larger is the tuning parameter, c∗ ⁠, the lower is the autocorrelation of the resulting posterior draws (though computation time increases because we are sampling with respect to more cluster locations). Good mixing is typically achieved with c∗ set equal to 2 or 3. We note that the number of clusters, M, may change in this step as clusters are added and deleted. 3. Sample DP concentration parameter, α: We use the algorithm of Escobar and West (1995) that formulates the posterior for α as a mixture of two Gamma distributions with the mixture component, η, drawn from a beta distribution. The algorithm is facilitated by the conditional independence of α from data, Y, given cluster assignments, s, (that implies number of clusters, M). 4. Draw fi as post-processing step from its posterior predictive distribution: use the Gaussian joint distribution for (fi,yi) to formulate the conditional posterior predictive Gaussian density, fi|yi∼NT(mi,Vi),(8) where mi=C(θsi∗)[Cτ(θsi∗)]−1yi and Vi=C(θsi∗)−C(θsi∗)[Cτ(θsi∗)]−1C(θsi∗) ⁠. 5. In the case of missing data values, we co-sample {fi} (rather than marginalizing over them) after the cluster locations in a Gibbs scan from, π(fi|yi,Θ∗,τɛ)=NT(φi−1ei,φi−1), where T×1, ei=fi′τɛyi and φi=τɛIT+C(θsi∗)−1 ⁠, where we cache the cholesky decompostion of C from sampling locations, Θ∗ ⁠, for faster computation of C−1 ⁠. Under co-sampling of the {fi} ⁠, we draw the model error precision, τɛ ⁠, in a Gibbs step (instead of in the Metropolis-Hastings scheme along with Θ∗ ⁠) using, π(τɛ|Y,{fi})=G(a1,b1), where shape hyperparameter, a1=0.5TN+a and rate hyperparameter, b1=0.5∑i=1N(yi−fi)′(yi−fi)+b and (a, b) are the prior hyperparameter settings. Lastly, we replace yi with fi under co-sampling of fi in (6) and (7) for sampling cluster locations, Θ* ⁠, and cluster memberships, s. A.2 Intrinsic Gaussian Markov Random Field Gibbs Scan Posterior sampling from the mixtures of iGMRFs employs a Gibbs scan over λ=({fi},{κm∗},s,α,τɛ) ⁠, from a set of conjugate full conditionals in the following steps: 1. Sample latent functions, {fij}i=1,…,N; j=1,…,T in a Gibbs steps from the full conditionals, π(fij|{fij}−j,τɛ,yij)=N(eijφij,φij−1), where eij=τɛyij+Qjjκsi∗f¯ij ⁠, with f¯ij=−1Qjj∑k:k∼jQjkfik and φij=τɛ+Qjjκsi∗ ⁠. 2. Sample locations, {κm∗} ⁠, in a Gibbs step from, π(κm*|{fij})=Ga(a2,b2), with shape parameter, a2=12nm(T−oQ)+a ⁠, where oQ = 2 is the rank-deficiency of the precision matrix, Q, indicating that that fi provides the equivalent of T−oQ degrees of freedom, rather than number of time points, T. The rate parameter, b2=12∑i:si=m(∑j=1TQjj[fij−f¯ij]2)+b ⁠, where the rate parameter is composed from the subset of latent functions, {fi}i:si=m ⁠, for those domains assigned to cluster m. 3. Sample cluster assignments under a similar Pólya urn representation as for the GP, only here the mixture posterior is conjugate, so, π(si=s|s−i,κ∗,α,τɛ,fi)∝{n−i,sn−1+α∏j=1TN(fij|f¯ij,[κsi∗Qjj]−1) if 1≤s≤M−αn−1+αd0,i if s=M−+1,(9) where d0,i=∫N(y|κ,…)G0(dκ)=baΓ(a2)[∏j=1JN(0|0,Qjj−1)]Γ(a)b2,ia2 ⁠, where a2 is as defined above and b2,i is term i in the sum that composes b2, defined above. Finally, we sample the DP concentration parameter, α, and the model noise precision, τɛ ⁠, in the same fashion as for the GP. Unlike the mixtures of GPs, we specify this model in a conjugate formulation that allows for fast sampling. B. ASSESSMENTS OF UNCERTAINTY IN CLUSTERING OF STATES IN CPS APPLICATION Figure 11 provides a 51 × 51 heat map of the empirical co-clustering probabilities used to summarize the posterior draws over the space of partitions for the CPS application of section 5. Relatively high and low values provide indication of higher concentration of the posterior distribution over the space of partitions. We observe relatively high and low values bounded away from 0.5. Figure 12 displays the MCMC trace plot of M, the total number of estimated clusters, that provides further insight in the uncertainty in the posterior distribution over the space of partitions. We observe both good mixing and a relative concentration around M = 2, which accords with our estimated MAP partition. Table 2 displays the distribution for values of M across the 1,000 post-burnin, thinned MCMC draws. Table 2. Distribution for Total Number of Estimated Clusters, M, Across the 1000 Post-Burnin MCMC Thinned Draws M . . 2 505 3 377 4 111 5 7 M . . 2 505 3 377 4 111 5 7 Open in new tab Table 2. Distribution for Total Number of Estimated Clusters, M, Across the 1000 Post-Burnin MCMC Thinned Draws M . . 2 505 3 377 4 111 5 7 M . . 2 505 3 377 4 111 5 7 Open in new tab Figure 12. Open in new tabDownload slide MCMC Trace Plot of the Number of Estimated Clusters, M, Over the Post-Burnin MCMC Draws. REFERENCES Banerjee S. , Wall M. M., Carlin B. P. ( 2003 ), “Frailty Modeling for Spatially Correlated Survival Data, with Application to Infant Mortality in Minnesota,” Biostatistics (Oxford) , 4 , 123 – 142 . Google Scholar Crossref Search ADS WorldCat Blackwell D. , MacQueen J. B. ( 1973 ), “Ferguson Distributions via Pólya Urn Schemes,” The Annals of Statistics , 1 , 353 – 355 . Google Scholar Crossref Search ADS WorldCat Dahl D. B. ( 2006 ), “Model-Based Clustering for Expression Data via a Dirichlet Process Mixture Model,” in ‘Bayesian Inference for Gene Expression and Proteomics’ , eds. Do K., Müller P., Vannucci M., pp. 201 – 215 , New York, NY : Cambridge University Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Escobar M. D. , West M. ( 1995 ), “Bayesian Density Estimation and Inference Using Mixtures,” Journal of American Statistical Association , 90 , 577 – 588 . Google Scholar Crossref Search ADS WorldCat Ferguson T. S. ( 1973 ), “A Bayesian Analysis of Some Nonparametric Problems,” The Annals of Statistics , 1 , 209 – 230 . Google Scholar Crossref Search ADS WorldCat Gelfand A. E. , Kottas A., Mac Eachern S. N. ( 2005 ), “Bayesian Nonparametric Spatial Modeling with Dirichlet Process Mixing,” Journal of the American Statistical Association , 100 , 1021 – 1035 . Google Scholar Crossref Search ADS WorldCat Lijoi A. , Mena R. H., Prünster I. ( 2007 ), “Controlling the Reinforcement in Bayesian Non-Parametric Mixture Models,” Journal of the Royal Statistical Society, Series B: Statistical Methodology , 69 , 715 – 740 . Google Scholar Crossref Search ADS WorldCat Neal R. ( 2000a ), “Slice Sampling,” Annals of Statistics , 31 , 705 – 767 . Google Scholar Crossref Search ADS WorldCat Neal R. M. ( 2000b ), “Markov Chain Sampling Methods for Dirichlet Process Mixture Models,” Journal of Computational and Graphical Statistics , 9 , 249 – 265 . Google Scholar OpenURL Placeholder Text WorldCat Pitman J. , Yor M. ( 1997 ), “The Two-Parameter Poisson-Dirichlet Distribution Derived from a Stable Subordinator,” The Annals of Probability , 25 , 855 – 900 . Google Scholar Crossref Search ADS WorldCat R Core Team ( 2014 ), R: A Language and Environment for Statistical Computing , R Foundation for Statistical Computing , Vienna, Austria . Available at http://www.R-project.org/. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rasmusen C. E. , Williams C. ( 2006 ), Gaussian Processes for Machine Learning , Cambridge, MA : The MIT Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rue H. , Held L. ( 2005 ), Gaussian Markov Random Fields: Theory and Applications , Boca Raton, FL : Chapman & Hall Ltd/CRC . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Savitsky T. D. , Vannucci M. ( 2010 ), “Spiked Dirichlet Process Priors for Generalized Gaussian Process Models,” Journal of Probability and Statistics , 2010 , 1 – 14 . Google Scholar Crossref Search ADS WorldCat Sethuraman J. ( 1994 ), “A Contructive Definition of Dirichlet Priors,” Stastica Sinica , 4 , 639 – 650 . Google Scholar OpenURL Placeholder Text WorldCat U.S. Bureau of Labor Statistics ( 2012 ), “BLS Spotlight on Statistics: The Recession of 20072009,” Available at https://www.bls.gov/spotlight/2012/recession/pdf/recession_bls_spotlight.pdf. (Accessed January 24, 2018). Vanhatalo J. , Riihimaki J., Hartikainen J., Jylanki P., Tolvanen V., Vehtari A. ( 2014 ), “Gpstuff,” Available at http://mloss.org/software/view/451/. (Accessed January 24, 2018). Wang C. , Neal R. M. ( 2013 ), “MCMC Methods for Gaussian Process Models Using Fast Approximations for the Likelihood,” Available at http://arxiv.org/abs/1305.2235v1. (Accessed January 24, 2018). Published by Oxford University Press on behalf of the American Association for Public Opinion Research 2018. This work is written by a US Government employee and is in the public domain in the US. Published by Oxford University Press on behalf of the American Association for Public Opinion Research 2018.
TI - Bayesian Nonparametric Functional Mixture Estimation for Time-Series Data, With Application to Estimation of State Employment Totals
JF - Journal of Survey Statistics and Methodology
DO - 10.1093/jssam/smy001
DA - 2019-03-01
UR - https://www.deepdyve.com/lp/oxford-university-press/bayesian-nonparametric-functional-mixture-estimation-for-time-series-H0LQv0WWj4
SP - 3
EP - 33
VL - 7
IS - 1
DP - DeepDyve
ER -