Spatio-Temporal Models for Big Multinomial Data using the Conditional Multivariate Logit-Beta Distribution
Jonathan R. Bradley, Christopher K. Wikle, Scott H. Holan
TL;DR
This paper presents a Bayesian spatio-temporal model for high-dimensional multinomial data using the conditional multivariate logit-beta distribution, enabling flexible covariance structures and efficient inference.
Contribution
It introduces the MN-STM, a novel Bayesian model that handles nonstationary, asymmetric covariances in space and time for multinomial data, with conjugate inference via the logit-beta distribution.
Findings
Model effectively captures complex spatio-temporal dependencies.
Demonstrates superior performance on simulated and real data.
Provides a scalable inference algorithm for high-dimensional multinomial data.
Abstract
We introduce a Bayesian approach for analyzing high-dimensional multinomial data that are referenced over space and time. In particular, the proportions associated with multinomial data are assumed to have a logit link to a latent spatio-temporal mixed effects model. This strategy allows for covariances that are nonstationarity in both space and time, asymmetric, and parsimonious. We also introduce the use of the conditional multivariate logit-beta distribution into the dependent multinomial data setting, which leads to conjugate full-conditional distributions for use in a collapsed Gibbs sampler. We refer to this model as the multinomial spatio-temporal mixed effects model (MN-STM). Additionally, we provide methodological developments including: the derivation of the associated full-conditional distributions, a relationship with a latent Gaussian process model, and the stability of the…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Chapter \thechapter
Spatio-Temporal Models for Big Multinomial Data using the Conditional Multivariate Logit-Beta Distribution
Jonathan R. Bradley111(to whom correspondence should be addressed) Department of Statistics, Florida State University, 117 N. Woodward Ave., Tallahassee, FL 32306-4330, [email protected], Christopher K. Wikle222Department of Statistics, University of Missouri, 146 Middlebush Hall, Columbia, MO 65211-6100, and Scott H. Holan2, 333U.S. Census Bureau, 4600 Solver Hill Road, Washington D. C. 20233-9100
Abstract
We introduce a Bayesian approach for analyzing high-dimensional multinomial data that are referenced over space and time. In particular, the proportions associated with multinomial data are assumed to have a logit link to a latent spatio-temporal mixed effects model. This strategy allows for covariances that are nonstationarity in both space and time, asymmetric, and parsimonious. We also introduce the use of the conditional multivariate logit-beta distribution into the dependent multinomial data setting, which leads to conjugate full-conditional distributions for use in a collapsed Gibbs sampler. We refer to this model as the multinomial spatio-temporal mixed effects model (MN-STM). Additionally, we provide methodological developments including: the derivation of the associated full-conditional distributions, a relationship with a latent Gaussian process model, and the stability of the non-stationary vector autoregressive model. We illustrate the MN-STM through simulations and through a demonstration with public-use Quarterly Workforce Indicators (QWI) data from the Longitudinal Employer Household Dynamics (LEHD) program of the U.S. Census Bureau.
Keywords: Bayesian hierarchical model; Big data; Pólya-Gamma; Markov chain Monte Carlo; Generalized Linear Mixed Model; Gibbs sampler.
1 Introduction
High-dimensional multinomial data referenced over space and time are ubiquitous among several disciplines. For example, machine learners are often interested in analyzing text corpora (or large sets of texts; Blei et al.,, 2003), demographers are interested in assessing populations among several categories across the U.S. (Hummer et al.,, 1995), ecologists often track the population of tree species over space and time (e.g., see Hobbs and Hooten,, 2015, Chp. 3), and epidemiologists monitor cancer rates by race over the U.S. (Chen et al.,, 2012). To better estimate proportions associated with these categories, one can leverage spatio-temporal dependence (e.g. see Cressie and Wikle,, 2011; Linderman et al.,, 2015, among others); however, spatio-temporal dependence can lead to difficult methodological and computational challenges. *Thus, the primary goal of this article is to develop a computationally feasible dynamic multinomial spatio-temporal model for high-dimensional datasets, with the purpose of estimating proportions. *
We are partially motivated by data from the Longitudinal Employer Household Dynamics (LEHD) program’s public-use Quarterly Workforce Indicators (QWIs). LEHD has become a leading authority on U.S. economics data (e.g., see Abowd et al., (2009), Abowd et al., (2013), and the references therein). We are motivated by the public-use QWI dataset (https://lehd.census.gov/data) partially because it is high-dimensional, has complex dependencies, and is an important tool for economists. Henceforth, we will use the phrase “complex dependencies” to indicate nonstationarity in both space and time, and asymmetry, which are well-known properties in the spatio-temporal literature (see, Cressie and Wikle, (2011) and Appendix A of the Supplementary Material of this article for a discussion of these terms). Additionally, many of the QWIs are not made available because some states do not sign the required Memorandum of Understanding (MOU) every year (Abowd et al.,, 2009, Sections 5.5.1 and 5.6). This suggests a need for small area estimation. Also, public-use LEHD currently does not release associated margins of error, and hence, there is an opportunity for to use a statistical model for uncertainty quantification.
LEHD provides data over all North American Industry Classification System (NAICS) sectors, all 3,145 U.S. counties, and over every yearly quarter. The public-use QWI dataset is very similar to high-dimensional panel data, or sometimes called longitudinal time series data (e.g., for recent references see Belloni et al.,, 2016; Gunawan et al.,, 2017a, and the references therein). However, a key difference is that models for panel data do not explicitly model spatial correlations, which are present in the LEHD dataset. We are particularly interested in estimating the proportion of individuals employed at the beginning of a quarter over each of the 20 NAICS sectors at each U.S. county. The size of this correlated dataset is quite large, at 2, 247, 586 observations.
It is important to be precise by what we mean by “big data.” For example, we are not considering the “ greater than ” problem, which is an important type of “big data” problem (Hastie et al.,, 2009). Instead, we are interested in the methodological difficulties involved with defining covariances when the sample size is large (e.g., millions of observations), which has become a common problem in the “big spatial data” literature (e.g., see Sun and Li,, 2012; Bradley et al.,, 2016; Heaton et al.,, 2018, for reviews).
There are several methods that one might adapt to analyze the high-dimensional LEHD dataset. In terms of dependent (spatial, spatio-temporal, or multivariate) multinomial data, there are many tools currently available. In particular, latent Dirichlet allocation (LDA; Blei et al.,, 2003) is a well-known model often used to solve the text corpora problem. LDA has been developed in the context of dependent (i.e., possibly in space and time) data, and is sometimes referred to as the correlated topic model (CTM; Blei and Lafferty,, 2005), which models dependencies by assuming that the logit of the proportions follow a multivariate normal distribution. The first implementation of the CTM (Blei and Lafferty,, 2005) used non-conjugate optimizations (e.g., see Andrieu et al.,, 2010; Murray et al.,, 2010). As a result, several recent papers have capitalized on the fact that the binomial distribution with a logit-link can be augmented using a Pólya-Gamma representation. This leads to easy to sample from conjugate full-conditional distributions (e.g., see Zhou et al.,, 2012; Chen et al.,, 2013; Polson and Scott,, 2011; Polson et al.,, 2013; Linderman et al.,, 2015; Glynn et al.,, 2018) for use in a Gibbs sampler. Additionally, this Pólya-gamma augmentation approach has been used to model dynamics (Blei and Lafferty,, 2006), through the incorporation of a vector autoregressive model (Blei and Lafferty,, 2006; Linderman et al.,, 2015).
The current state-of-the-art in the CTM literature requires one to augment the multinomial random vector with Póly-gamma random variables. These Póly-gamma random random variables are simulated using an algorithm with iterative calculations and an accept/reject step (Polson and Scott,, 2011), which is not computationally advantageous when doing repeated simulations within a Gibbs sampler. Thus, the first contribution of this article is to propose the use of the conditional multivariate logit-beta distribution (MLB) from Bradley et al., (2018a) to model dependent multinomial data. The MLB distribution gives conjugate full-conditional distributions that are straightforward to simulate from, and allows one to avoid computationally less efficient data augmentation. We provide technical a result that develops the relationship between the MLB distribution and other multivariate distributions. Most notably, one can view our hierarchical model as a special type of latent Gaussian process model. This MLB distribution has been recently introduced by Bradley et al., (2018a) and has been used to model Bernoulli, binomial, and negative binomial data in the multivariate/spatial settings. However, the MLB distribution has not yet been used in the multinomial spatio-temporal setting.
The computational advantages of the MLB distribution for modeling multinomial data are especially powerful when considering the current state of the literature for non-Gaussian dependent data (e.g., see Diggle et al.,, 1998, for standard early references.). There are many Bayesian methods available to analyze non-Gaussian data. Quite often these algorithms require difficult to tune Metropolis-Hasting steps nested within a Gibbs sampler (Shaby and Wells,, 2011; Roberts and Tweedie,, 1996; Garthwaite et al.,, 2010). There are methods that aid in tuning a Markov chain Monte Carlo (MCMC) algorithm, such as pseudo-marginal MCMC (Andrieu and Roberts,, 2009); however, these approaches have not been extended to the multinomial spatio-temporal setting. Tuning steps are completely avoided by our proposed approach, and by the method in Linderman et al., (2015). However, the method in Linderman et al., (2015) has not been developed for high-dimensional settings (e.g., it does not incorporate dimension reduction). Other algorithms such as Hamiltonian Markov chain Monte Carlo (e.g., see Neal,, 2003; Dang et al.,, 2017), particle MCMC (Gunawan et al.,, 2017b), and splice sampling (Murray et al.,, 2010) are also known to be difficult to implement when using high-dimensional datasets (e.g., see discussions in Rue et al.,, 2009; Bradley et al.,, 2018b).
There are, of course, other perspectives for non-Gaussian dependent data that exist outside the realm of a fully Bayesian analysis, which are not considered in this manuscript. For example, in the spatial/non-Gaussian setting, empirical Bayesian methods are available. Sengupta and Cressie, (2013b), Sengupta and Cressie, (2013a), and Sengupta et al., (2016) employ an empirical Bayesian algorithm that involves a Newton-Raphson algorithm, nested within an expectation maximization algorithm, which is then nested within a Markov chain Monte Carlo algorithm. Sengupta and Cressie, (2013b) provide empirical results that demonstrates that their empirical Bayesian analysis can lead to improvements in prediction over a specific fully Bayesian model and is computationally advantageous. However, empirical Bayesian approaches do not explicitly model the variability introduced by estimating process model parameters. Approximate Bayesian methods are also available for spatial and spatio-temporal settings including integrated nested Laplace approximations (INLA) (e.g., see Rue et al.,, 2009; Bakka et al.,, 2018, and the references therein) and variational Bayes methods (Ren et al.,, 2011; Nathoo et al.,, 2014). However, these models have not been extended to the high-dimensional multinomial spatio-temporal setting that we are considering.
In addition to our new use of the MLB distribution, we also provide an efficient parameterization of the spatio-temporal covariance. In particular, we use two commonly used techniques to parameterize spatial/spatio-temporal covariances. The first technique, is to allow different random variables to share the same random effect. These techniques have been used in the spatial (e.g., see Banerjee et al.,, 2008; Cressie and Johannesson,, 2008; Lindgren et al.,, 2011; Hughes and Haran,, 2013; Nychka et al.,, 2014; Sengupta and Cressie,, 2013b, a; Sengupta et al.,, 2016; Bradley et al.,, 2018a, among others), multivariate-spatial (Royle et al.,, 1999; Finley et al.,, 2009, 2010), and spatio-temporal (Waller et al.,, 1997; Wikle et al.,, 2001; Bradley et al.,, 2015b) settings. The second commonly used technique is to partition the joint likelihood into a product of more manageable conditional likelihoods. For example, see the early papers on conditional autoregressive (CAR) models by Besag, (1974), Besag, (1986), and Besag et al., (1991) in the spatial setting; Royle and Berliner, (1999), Mardia, (1988), and Billheimer et al., (1997) for the multivariate spatial setting; and Wikle and Cressie, (1999), Cressie et al., (2010), Katzfuss and Cressie, (2011), and Katzfuss and Cressie, (2012) in the spatio-temporal setting.
Both techniques are used to define the multivariate spatio-temporal mixed effects (MSTM) introduced in Bradley et al., (2015a). In particular, to incorporate multivariate-spatial dependencies the MSTM lets the different categories and spatial regions share the same random effects at time (i.e., the first technique). Then, to incorporate temporal dependence the MSTM uses an order one vector autoregressive model VAR(1) (i.e., the second technique). The MSTM was first introduced for Gaussian data in Bradley et al., (2015a), and later adapted to Poisson data (Bradley et al.,, 2018b). Our second contribution is to adapt the MSTM to the multinomial data setting. The main methodological development needed here is to extend prior distributions of precision parameters (i.e., inverse of the covariance matrix) when using random effects distributed according to the MLB distribution. Previous MSTMs specify the model based precision parameter of the random effect to be “close” (in terms of the Frobenious norm) to the precision parameter implied by a CAR model. We use the same strategy in this article, and call our proposed model the multinomial spatio-temporal mixed effects model (MN-STM).
An important motivator of the MSTM is that it allows for dimension reduction, which aids in analyzing high-dimensional data (such as the LEHD dataset). Thus, our third contribution is to use a reduced rank model for multinomial spatio-temporal data, which has not yet been proposed in the CTM literature. Dimension reduction methods are often used when there is a computational bottleneck when computing a high-dimensional likelihood. To address this issue, latent high-dimensional random vectors are replaced by low-dimensional random vectors; in multivariate analysis, this is similar to principal component analysis (e.g., see Jolliffe,, 2002; Cox,, 2005, among others). Dimension reduction in the spatial/spatio-temporal setting has a mature literature associated with it. For example, see Cressie and Johannesson, (2006), Cressie and Johannesson, (2008), Shi and Cressie, (2007), Banerjee et al., (2008), Kang and Cressie, (2011), Lindgren et al., (2011), Nychka et al., (2014) for the univariate-spatial setting; and Wikle and Cressie, (1999), Cressie et al., (2010), Katzfuss and Cressie, (2011), and Katzfuss and Cressie, (2012) in the spatio-temporal setting. Also see, Wikle, (2010) for a discussion on spatial and spatio-temporal dimension reduction modeling.
In addition to dimension reduction, the MSTM also adapts aspects of the model suggested by Hughes and Haran, (2013) from the spatial only setting to the multivariate-spatio-temporal setting. Specifically, the propagator matrix associated with a VAR(1) model is specified so that random effects are not confounded (Bradley et al.,, 2015a, 2018b). This propagator matrix is referred to as the Moran’s I (MI) propagator matrix because of a connection to the MI statistic from Moran, (1950). This motivation is similar to the specification of the MI basis functions used in Griffith, (2000), Griffith, (2002), Griffith, (2004), Tiefelsdorf and Griffith, (2007), Hughes and Haran, (2013), Porter et al., (2013), and Burden et al., (2015), which we review in this manuscript. The dynamic properties associated with the MI propagator matrix still need development. In particular, when introducing a VAR(1) model for forecasting, it is especially important to verify that the Wold representation of the VAR(1) model exists (e.g., see Anderson,, 1971, for a standard reference). However, Bradley et al., (2015a) and Bradley et al., (2018b) did not investigate this. Thus, our fourth methodological contribution is to provide the conditions in which the Wold representation of the VAR(1) model (with MI propagator matrices) exists.
The remainder of this article is organized as follows. In Section 2, we introduce the use of the MLB distribution to model multinomial data with fixed and random effects. We end Section 2 with an illustrative example of logistic regression with MLB random effects, which is used to motivate the MN-STM to model data with more complex features (e.g., LEHD). In Section 3, we define the MN-STM and present the necessary technical development. We also provide a small technical result showing a type of equivalence between the use of the MLB distribution and the Gaussian distribution. This is done in an effort to “de-mystify” our use of new distribution theory. Empirical results are discussed in Section 4. This includes a simulation study that illustrates the computational performance of the MN-STM relative to a latent Gaussian process model and the the Poisson multivariate spatio-temporal mixed effects model (PMSTM) from Bradley et al., (2018b). Additionally, an example analysis using a dataset obtained from LEHD is offered. Finally, Section 5 contains a discussion. For convenience of exposition, some statements and proofs of technical results are given in the appendices.
2 Logistic Regression with Latent Multivariate Logit-Beta Random Effects
Consider multinomial data that are recorded at small areas and at equally spaced discrete time-points, where we assume that is empty for . Let be the multinomial data vector with categories, observed at time and areal unit , where each is integer-valued from zero to . For example, might consist of counts of the number of people employed in different NAICS sectors at county and yearly-quarter . The goal of our analysis is to estimate the probability of each category for every and . Let represent the counties at time where a multinomial vector is observed, and suppose there are areal units in and let .
2.1 Multinomial Data
We assume that follows a multinomial distribution. That is,
[TABLE]
where it is assumed that and for every and . Now, write the pmf of the multinomial distribution as the product of different binomial probability mass functions (pmfs) as follows (e.g., see Linderman et al.,, 2015, among others):
[TABLE]
where , and . Let , so that Equation (1) can be re-expressed as
[TABLE]
Hence, the pmf for the -dimensional vector is given by,
[TABLE]
where the -dimensional vector , the -dimensional vector , and the symbol “” can be read as “proportional to as a function of .”
2.2 The Multivariate Logit-Beta Distribution
The derivation of the MLB distribution starts with the univariate logit-beta random variable, which is defined as
[TABLE]
where is a beta random variable with first shape parameter and second shape parameter . The shape parameter is parameterized so that it is strictly larger than . This choice of parameterization will simplify subsequent expressions. The transformation in (4) leads to the following probability density function (pdf) for ,
[TABLE]
where recall that is used to denote a generic pdf. Then, the MLB distribution is found through the following transformation:
[TABLE]
where is an -dimensional real-valued vector, V is an real-valued invertible matrix, , and the ’s are independent logit-beta random variables with shapes and , respectively. Here, is a vector-valued location parameter and V is a matrix-valued precision parameter. Again, straightforward change-of-variables of the transformation in (6) yields the following expression for the pdf of q (Bradley et al.,, 2018a):
[TABLE]
where is a -dimensional vector of ones, , and . The pdf in (7) has a similar functional form as the likelihood associated with a multinomial probability mass function with a logit-link in (3), and is the aforementioned MLB distribution. This relationship leads to conjugacy between the MLB distribution and the multinomial distribution.
The conditional MLB distribution also has a similar form as (7). Denote the partioning of the vectors and , such that and are -dimensional, and and are -dimensional. Also, partition into the matrix H and the matrix B. Then the conditional MLB distribution is given by
[TABLE]
where d is a fixed real-valued -dimensional vector. Define the -dimensional vector . For ease of notation, we will let
[TABLE]
Here, c can be equal to any -dimensional real-valued vector . Specifically, re-parameterize so that and . This will be useful later on, as the full-conditional distributions associated with our model will have this re-parameterized form. Furthermore, we drop and and include c in our notation, and signify this notational convenience with the subscript “MLB” in , and refer to (9) as the conditional MLB distribution. We refer to c as the location parameter.
It is difficult to directly simulate from a conditional MLB. However, in a Bayesian context one can augment the likelihood in a manner so that updates within a Gibbs sampler can be computed using marginal distributions from a MLB. In particular, we will use the following result.
Proposition 1: Let , where is full column rank, , , , and for . Assume a re-parameterized value (see discussion below (9)) of , and the improper prior , where is -dimensional. Also let be the orthonormal basis for the null space of , , , be an identity matrix, and let . Then,
[TABLE]
*where the integrand in (10) is proportional to for defined in (7). Furthermore, the affine transformation, *
[TABLE]
*is a draw from the density in (10).
Proof: See Appendix A.
Proposition 1 will allow us to augment the likelihood and update multinomial parameters using a collapsed Gibbs sampler (Liu,, 1994). That is, when a full-conditional distribution is proportional to a we will collapse across the mean, which will be assumed to have an improper prior. We provide an example of this in a simplified setting in Section 2.3.
2.3 Example: Logistic Regression with Latent Multivariate Logit-Beta Random Effects
We now provide an illustration of Proposition 1 for implementing a logistic regression model with latent multivariate logit-beta random effects. Consider and the following hierarchical model:
[TABLE]
where , the is a -dimensional vector of covariate effects, is a -dimensional vector of random effects, is a -dimensional vector of random effects, , , and . Let be an identity matrix, be a -dimensional vector of zeros, let the -dimensional vector consist of strictly positive elements, and the -dimensional vector have elements strictly larger than the corresponding elements of . For illustration, let X be an matrix of orthonormal covariates and let be an orthonormal “design matrix” (Hodges,, 2013). The prior distributions in (32) depends on the sample size (e.g., the precision parameter for the prior of is -dimensional). Thus, we introduce and choose it to be small to mitigate changes in the prior as more samples are collected. The need for both and is discussed at the end of Section 2.3.
This hierarchical model might be implemented with the Gibbs sampler, which is outlined in Pseudo-Code (1). To implement Pseudo-Code (1) we first find the full conditional distribution for as follows,
[TABLE]
where and . Similarly, the full-conditional distributions for and are given by,
[TABLE]
As discussed in Section 2.2, it is difficult to directly sample from these full-conditional distributions; hence we offer a data augmentation scheme that leads to full-conditional distributions that are easy to sample from. To do this, we augment the likelihood in (32) with the -dimensional random vector , the -dimensional random vector , and the -dimensional random vector . Now, consider the following hierarchical model:
[TABLE]
where is the orthonormal basis of the orthogonal matrix , where is , is , and is . Similarly, is the orthonormal basis of the orthogonal matrix , where is , is , and is . Likewise, is the orthonormal basis of the orthogonal matrix , where is , is , and is .
In (86), the conditional distribution of and given , , and is proportional to the likelihood in (32). That is, in (86) replace with , with , and with to obtain a likelihood proportional to (32). Thus, a strategy to sample from a likelihood proportional to (32), would be to sample from . This can be easily implemented using the following collapsed Gibbs sampler (Liu,, 1994) as outlined in Pseudo-Code (2).
In Pseudo-Code 2, we have collapsed across the event in Step 2, in Step 3, and in Step 4. To do this, we first find and then marginalize across . That is,
[TABLE]
It follows from Proposition 1 that is a sample from , where w is distributed as . Similar algebra shows that is a sample from , where w is distributed as . Likewise, is a sample from , where w is distributed as . Figure (1) shows an example of and the corresponding predicted values from data simulated with this choice of (details behind this small simulation are given in Appendix B). The estimated values appear to be trending the true values quite well in this small illustration.
Let . Then, it is possible to have a shape parameter equal to zero, since it is possible to observe a zero count. Similarly, if and , it is possible to have the two shape parameters equal to each other, because it is possible to observe a . Thus, the careful attention required to introduce strictly positive , , and allows us to use Proposition 1 to directly sample in a collapsed Gibbs sampler.
This example model illustrates the need for , , , and Proposition 1. However, a consequence of this is that there is a considerable amount of bookkeeping needed. Despite the amount of attention required to implement such a collapsed Gibbs sampler, this approach allows one to directly sample from the full-conditional distribution in straightforward manner; and hence, avoids tuning proposal densities in a Metropolis-Hastings algorithm. We demonstrates the computational gains of this type of approach in Section 4.1.
3 Spatio-Temporal Logistic Regression with Latent Multivariate Logit-Beta Random Effects
The incorporation of low-dimensional MLB random is a significant contribution to the CTM literature. However, modern datasets often exhibit more features than what the illustrative model in Section 2.3 can capture. In particular, many high-dimensional datasets exhibit dependent variation over time (i.e., dynamics). As such, additional model details are required to match this feature often found in modern datasets, and hence, in Section 3 we offer a dynamic modeling version of the illustrative model presented in Section 2.3. We refer to this dynamic model as the MN-STM. A complete summary of the MN-STM, and additional discussion on basic properties of the MN-STM are presented in Appendix A of the Supplementary Material. The full-conditional distributions associated with the MN-STM are derived in Appendix B of the Supplementary Material.
3.1 The Process Model
Similar to Section 2.3, we make the following assumption on :
[TABLE]
where the -dimensional vector of covariates are assumed known, and is an associated unknown -dimensional vector. The -dimensional real-valued vector is known, and is an -dimensional random vector with elements that are possibly correlated. The elements of are assumed to be independent and identically distributed.
The prior for the dimensional vector is assumed to be with location-parameter-zero; has unknown shape parameters, , , , , and the elements of the -dimensional vector are positive, where ; and precision parameter , where the matrix , the block diagonal matrix , and is a identity matrix. The -dimensional random vector is assumed to be with location-parameter and unknown shape parameters , , and , where , the -dimensional vectors stack to produce , the precision parameter , and is unknown. For a given , define the -dimensional vector follow a with zero location parameter, precision parameter , and constant (across and ) shape parameters and .
The specification of is extremely important. If is chosen poorly, then the fixed effects and the random effects may be confounded (e.g., see Wilson and Reich,, 2014, and the references therein). Consider the extreme case where . Under this specification there is no way to disentangle the fixed effect and the random effect when estimating the expected values of . One solution to this problem is to specify so that it belongs to the orthogonal column space of (see, Griffith,, 2000, 2002, 2004; Tiefelsdorf and Griffith,, 2007). An approach used for areal data is to consider a spatially weighted version of the orthogonal complement. This is referred to as the so-called “Moran’s I” basis functions (see, Hughes and Haran,, 2013; Porter et al.,, 2013; Burden et al.,, 2015; Bradley et al.,, 2015a). This spatially weighted orthogonal complement is the classical Moran’s I operator (see, Moran,, 1950; Hughes and Haran,, 2013). Specifically, let the matrix , where “P” represents “prediction locations” and delineates between which only stacks the covariates over observed regions. The Moran’s I operator is written as,
[TABLE]
where is an identity matrix, and A is the adjacency matrix associated with the edges formed by . One choice for the adjacency matrix is to set the -th element of the matrix A equal to one if the areal unit is a neighbor of (Cressie,, 1993; Banerjee et al.,, 2015). Let be the spectral decomposition of the matrix . Then set equal to the first columns of .
For a given and the model in Equation (87) is a type of spatial mixed effects model (e.g., see Cressie and Johannesson,, 2008, among others). The random vector is the “shared random effect” discussed in the Introduction, meaning the expressions of and using Equation (87) both contain the same random vector for and . This shared random effect induces cross correlations among different categories (i.e., and ) and different spatial regions (i.e., and ), since
[TABLE]
which is not equal to zero. As discussed in the Introduction, there are several examples of this technique in the spatial, multivariate spatial, spatio-temporal, and multivariate spatio-temporal literature.
3.2 Spatio-Temporal Dynamics
Both Blei and Lafferty, (2006) and Linderman et al., (2015) consider the use of a VAR(1) model to incorporate dynamics within multinomial data. In this article, we consider incorporating a specific type of VAR(1) model used in Bradley et al., (2015a) and Bradley et al., (2018b).
The VAR(1) assumption for the -dimensional non-spatially referenced random vector is given by (e.g., see Cressie and Wikle,, 2011, Chap. 7),
[TABLE]
where is a known propagator matrix (see discussion below), and is an -dimensional random vector. We use the so-called Moran’s I propagator matrix introduced in Bradley et al., (2015a). The derivation of this matrix starts by substituting the VAR(1) representation into the mixed effects model in (87), at all prediction locations, to obtain
[TABLE]
Then, straightforward algebra leads to,
[TABLE]
where the matrix and the -dimensional random vector . Similar to the derivation of , we let equal the eigenvectors of , where in general is a real-valued “weight” matrix. In Section 5, we set .
When introducing a VAR(1) model for forecasting, it is especially important to verify that the Wold representation of the VAR(1) exists (e.g., see Anderson,, 1971, for a standard reference). This implies that the VAR(1) model is stable. The word “stable” means that the variance of the random effects does not drift to infinity as increases. Our primary use of the VAR(1) model is to incorporate time-dependence to aid in predicting over the range times that were observed (i.e., smoothing), and not forecasting. However, it may be of interest to use our VAR(1) model for forecasting in future research. Thus, we provide a result that gives the necessary conditions for the MI propagator to be used to define a stable process (e.g., see Lutkepohl,, 2005, pg. 688).
*Proposition 2: Consider the VAR(1) model in (89). Denote the first eigenvectors of with and let be a generic real-valued matrix. For , let . Also, for let for some . Then, the Wold representation of the VAR(1) model in (89) converges in to a stable sequence of -dimensional random vectors.
Proof: See Appendix A.
In Section 5, we set because we are only interested in smoothing. However, if forecasting is of interest, one might place a prior distribution on and set .
3.3 Marginal Covariance of the Random Effects
For computational reasons, we do not include all random effects (i.e., ) and confounded random effects in (87). Thus, we choose so that the spatially dependent term in (87), namely , has a precision (i.e., inverse of the covariance matrix) similar to the precision from an intrinsic conditional autoregressive (ICAR) model, which is given by for . It is well-known that ICAR models induces covariances that are functions of neighborhoods and are not functions of distances between locations; thus, the ICAR model enforces nonstationary spatial dependence (Besag,, 1974).
We let , where
[TABLE]
For a square real-valued matrix C let be the generalized inverse of C and let the Frobenius norm . In (91), we minimize the Frobenius norm across the space of all real-valued matrices. The general strategy in (91) is that same as in Bradley et al., (2015a) and Bradley et al., (2018b); however, the term is different from the associated covariance terms in Bradley et al., (2015a) and Bradley et al., (2018b) because we use the MLB distribution.
Proposition 3: Let the model in (87) hold, and let P be a generic positive semi-definite matrix. Then the values of and that minimize,
[TABLE]
are given by:
[TABLE]
*where , is the best positive approximate Higham, (1988) of the matrix C, is the spectral decomposition of , and is the trigamma function (Johnson,, 1949).
Proof: See Appendix A.
We set P in (92) to . One could also choose to specify the joint precision of to be close to a spatio-temporal autoregressive model. However, the act of stacking the basis functions over quickly leads to memory issues. Consequently, we specify marginally over each using Proposition 3. Additionally, Proposition 3 shows that the shape parameters define the variances of the random effects , and consequently, we apply gamma priors to every shape parameter. Thus we include shape parameter priors for and (see Appendix A of the Supplementary Material for more details).
3.4 A Latent Gaussian Process Model Representation of the MN-STM
The latent Gaussian process model is a well established method used to analyze spatial and spatio-temporal data (e.g., see Gelfand and Schliep,, 2016, among others), and hence, one might be hesitant to adopt a new distributional framework. Motivated by this, we show that the latent MLB model can be written as a latent Gaussian process model. Our strategy is to augment the MLB likelihood using Póly-gamma random variables similar to Polson et al., (2013). The main difference between our result and Polson et al., (2013) is that Polson et al., (2013) augments the binomial pmf, while we augment the MLB distribution.
Proposition 4: Let , , and , where is a Póly-gamma random variable with density . The definition of is given in Appendix A. Denote the joint pdf of the independent Póly-gamma random variables with . Then,
[TABLE]
*where is a strictly postive real-valued function of and .
Proof: See Appendix A.
Proposition 4 suggests that the MLB distribution can be seen as an infinite mixture of normal densities (i.e., using the Reimann sum expression of the integral). This implies that the MN-STM can be augmented by random variables so that it can be re-expressed as a latent Gaussian process model (see Appendix C of the Supplementary Material for this expression).
4 Empirical Results
4.1 A Simulation Study
Both Bradley et al., (2015a) and Bradley et al., (2018b) conduct an empirical simulation analysis to understand the out-sample performance of their model if they were able to obtain independent replicates of QWIs. Specifically, the unknown parameters in the data model are replaced by public-use QWIs. Simulated values are then generated from this empirical data model, and the simulated data are used for prediction. The main motivation of this approach is that it yields simulated values that are similar to what might be seen in practice. We adopt this strategy to build our simulation model.
Let be the beginning of quarter employment QWI, at quarter , for industry , and Minnesota county . We simulate data, denoted by , from a multinomial distribution with probability for and . Let consist of counties in Minnesota (MN) that have available QWIs. Here, denotes the information industry, represents the professional, scientific and technical services industry, and represents the finance industry. We add 1 in the numerator and 3 in the denominator of the probability of success so that the proportions are greater than zero and sum to one. We randomly select 65 of the areal units in to be “observed.” The covariates , where is the indicator function. These covariates were the same used in the simulation study in Bradley et al., (2018b). We set equal to approximately the top of the available basis functions. In this case, we set .
In Figure (2), we present sample proportions using the QWIs, the simulated data, and the predictor (see Section 3.5). The predictions are fairly accurate (i.e., the predictions are close to the sample proportions). The high quality of these predictions are further supported by Figure (3), where we display a scatterplot of the true proportions versus the estimated proportions.
We repeated this example over 50 independent replications of the set , and compare to the performance of the estimates using the MN-STM, PMSTM, and the latent Gaussian process (LGP) model considered in Bradley et al., (2018b). In Figure (4) we plot boxplots of the median relative absolute error,
[TABLE]
where “median” represents the median function over , , and (including , , and for unobserved ), “abs” denotes the absolute value function, and represents a generic estimate of . Values close to one show that the estimated value produces estimates that are comparable to the variability of the data at observed locations, and values less (greater) than one show an increased (decreased) performance relative to the variability of the data at observed locations. In Figure (4), we see that both the latent MLB (LMLB) model introduced in Section 2.3, and the P-MSTM performs reasonably well in terms of median relative absolute error with values close 1, while the LGP performs considerably worse than the P-MSTM. The MN-STM clearly outperforms the three competitors with median absolute errors close to 0.4. Both the PMSTM and the MN-STM appear to outperform the MLB, which does not explicitly model the dynamics through a VAR(1) model.
The computational performance of the MN-STM is also important to investigate. We use the effective sample size (ESS) to assess the performance of the MCMC chain, which is computed as the length of the MCMC () times the ratio of the within chain variance and the between chain variance. Small (large) values of ESS suggest that the MCMC chain has positive (negative) correlations between values in the chain. This has allowed many to use ESS as a measure of the efficiency of the MCMC, since ESS close to or larger than the length of the MCMC implies an efficient MCMC (standard references include Kass et al., (2016), Liu, (2008), Robert and Casella, (2013), and Gong and Flegal, (2016) for component-wise ESS, and Vats et al., (2016) for a multivariate ESS). In Figure (5), we show a boxplot of ESS ( over , , and ), for a single run of the MCMC (post burnin, ) and for a single generation of . The component-wise ESS of MN-STM and P-MSTM are comparable with the MN-STM outperforming the P-MSTM (the median ESS values are 899 and 430, respectively), and indicates computationally efficient Markov chains. The ESS of the LGP model is strikingly smaller than 1,000 suggesting that the Metropolis-within-Gibbs implementation of this LGP is extremely inefficient (the median ESS is approximately 19).
4.2 Application: Predicting the Probability of Employment in a NAICS Sector
We now compute predictions of the probability of employment in a NAICS Sector. The main purpose of this demonstration is to show that it is possible to jointly analyze such a high-dimensional dataset. We observe data over all NAICS sectors, U.S. counties, and quarters, using the high-dimensional QWI dataset of size 2, 247, 586 (see Figure (2)). We use the same covariates from Bradley et al., (2018b) and fit the Gibbs sampler based on the full-conditionals given in Appendix B of the Supplementary Material. Specifically, we let , where is the 2010 decennial Census value of the population of county and is the indicator function. We consider , where the first 100 eigenvalues of the Moran’s I operator makes up 99.5 of the total variability of the Moran’s I operator. The Gibbs sampler ran for 10,000 iterations, and had a burn-in of 1,000. Trace plots of the sample chains were used to check for convergence of the Gibbs sampler, and no lack of convergence was detected. The latent Gaussian process model did not mix well enough to provide any meaningful comparisons, and the model in Section 2.3 lead to storage errors.
In Figure (6), we (partially) plot the sample proportion of individuals employed at the beginning of the quarter, and the corresponding predictions and posterior standard deviation. In Figure (6), we see that the predictions roughly track the patterns of the LEHD (QWI) proportions. We are also able to produce measures of uncertainty (i.e., posterior standard deviations), which are also presented in Figure (6). These results are consistent across different subsets of the predictions.
5 Discussion
We have introduced methodology for jointly modeling dependent multinomial spatio-temporal data within the Bayesian framework. We refer to our model as the multinomial spatio-temporal mixed effects model (or MN-STM), which provides an avenue to model complex dependencies including nonstationarity, non-separability, and asymmetries. The proposed model is extremely parsimonious as we incorporate the Moran’s I propagator and introduce a marginal multivariate-spatial precision chosen to be close to an intrinsically conditional autoregressive model. We show that this VAR(1) model is stable, which is an important contribution to the Moran’s I propagator literature. Furthermore, the MN-STM is a computationally efficient, since the random effects are projected onto a reduced dimensional space.
A critical contribution of the MN-STM is the use of the MLB distribution (Bradley et al.,, 2018a) offers an exciting model to consider in the correlated topics model literature (Blei and Lafferty,, 2005, 2006; Linderman et al.,, 2015). The MLB distribution is conjugate, which allows one to avoid expensive data augmentation schemes. In particular, the introduction of a Póly-gamma random variable is not necessary to obtain conjugacy. We show that the MLB can be expressed as an infinite mixture of multivariate normal distributions, and hence can be expressed as a more conventional latent Gaussian process model.
Simulation studies showed that the MN-STM based on latent MLB random vectors performs extremely well. Specifically, the simulated examples all indicated very small out-of-sample error of the MN-STM model, and outperforms the P-MSTM and an latent Gaussian process model. Additionally, the computational performance of the MN-STM was shown to be efficient in terms of effective sample size. Furthermore, the MN-STM performed quite well even though the simulated data were not generated from the MN-STM, which suggests that our model is robust to departures from model assumptions. We illustrated that the MN-STM is scalable to the high-dimensional settings by analyzing a big count-valued dataset (of 2, 247, 586 observations) consisting of public-use QWIs made available by the U.S. Census Bureau’s LEHD program.
Acknowledgments
This research was partially supported by the U.S. National Science Foundation (NSF) and the U.S. Census Bureau under NSF grant SES-1132031, funded through the NSF-Census Research Network (NCRN) program. This article is released to inform interested parties of ongoing research and to encourage discussion. The views expressed on statistical issues are those of the authors and not necessarily those of the NSF or U.S. Census Bureau.
Appendix A: Proofs
Proof of Proposition 1: The distribution of is equal to , where recall we have reparameterized (see the discussion below Equation (9)). Thus,
[TABLE]
Integrating out we obtain,
[TABLE]
Thus, is the marginal random vector associated with for defined in (7). Thus, we are left to show that is a sample from this marginal distribution.
Denote the QR decomposition of , where the matrix Q satisfies and R is a upper triangular matrix. Now recall the definition of the matrix , which satisfies and . Then can be written as
[TABLE]
It follows that
[TABLE]
where the last equality in the above can be verified by substituting into . From (6), distributed according to can be written as
[TABLE]
where the -dimensional random vector w is distributed according to where is defined in (7). Multiplying both sides of (A.6) by we have
[TABLE]
and hence the distribution associated with is the marginal distribution associated with
as desired.
Proof of Proposition 2: Without loss of generality let . Then iterating the VAR(1) expression gives:
[TABLE]
where
[TABLE]
For the Wold representation to hold, we have from Fuller, (1976, pg., 29 31) that we need to show two items:
[TABLE]
To show (A.8), we have
[TABLE]
where for every , since is orthonormal for every . Now, upon taking the limit at goes to infinity,
[TABLE]
To show (A.9), first recall that the mean of a logit-beta random variable is , where is the digamma function. Also, the variance is given by , where is the trigamma function (Johnson,, 1949). From (11),
[TABLE]
where the consists i.i.d. logit-beta random variables with shape and . Thus, , where . Thus,
[TABLE]
provided that is invertible.
Proof of Proposition 3: Augment with an -dimensional random vector as done in Proposition 1, where the matrix be the orthogonal complement of . Then, from (A.6),
[TABLE]
where consists of i.i.d. logit-beta random variables with the -th element of is denoted with . It follows that . This gives us,
[TABLE]
Hence,
[TABLE]
Substituting (A.11) into (A.10), and a few lines of algebra, leads to
[TABLE]
Substituting (A.12) into (92), we obtain
[TABLE]
Let so that minimizing (A.13) is that same as minimizing,
[TABLE]
It follows from Proposition 1 in Bradley et al., (2015a) that the value of that minimizes (A.14) is
[TABLE]
Then, it follows that minimizes the Frobenius norm in (92). Finally, it is well known that (Johnson,, 1949)
[TABLE]
which completes the proof.
Proof of Proposition 4: Let . Then the Pólya-Gamma random variable is defined to be,
[TABLE]
has a density , which does not have a closed-form expression. We make use of the follow integral identity (e.g., see Polson et al.,, 2013):
[TABLE]
Let . Denote the normalizing constant of with , which is known to be finite (Bradley et al.,, 2018a). Then,
[TABLE]
which after a few lines of algebra becomes,
[TABLE]
where,
[TABLE]
This completes the result. As a small side-note, notice that implies that
[TABLE]
where the expectation is taken with respect to .
Appendix B: Specifications of the Small Simulated Example in Section 2.3
To produce Figure (1), we set = 1, , , and . We define to be the orthogonalization of and set equal to the first columns of the orthogonal complement of X. Set , the elements of to be values generated from a normal distribution with mean zero and variance 1, and generate the elements of to be from a normal distribution with mean zero and variance 0.33.
We run the MCMC for iterations, and trace plots do not indicate a lack of convergence. The value of , , , , , and . In our experience, the performance of Pseudo-Code 2 improves as the elements of get close to . Thus, our specification of might be interpreted as an empirical Bayes specification of Pseudo-Code 2, where we substitute an estimate of .
Chapter \thechapter
Supplementary Appendix: Spatio-Temporal Models for Big Multinomial Data using the Conditional Multivariate Logit-Beta Distribution
Jonathan R. Bradley444(to whom correspondence should be addressed) Department of Statistics, Florida State University, 117 N. Woodward Ave., Tallahassee, FL 32306-4330, [email protected], Christopher K. Wikle555Department of Statistics, University of Missouri, 146 Middlebush Hall, Columbia, MO 65211-6100, and Scott H. Holan2, 666U.S. Census Bureau, 4600 Solver Hill Road, Washington D. C. 20233-9100
Introduction
In this Supplementary Appendix, a summary and further description of the MN-STM (Supplementary Appendix A), the derivation of the Gibbs sampler for the MN-STM (Supplementary Appendix B), and an alternative latent Gaussian process expression of the MN-STM (Supplementary Appendix C).
Appendix A: Summary of the MN-STM
We choose to write the model for using the “data model,” “process model,” and “parameter model” notation that is commonly used in the spatio-temporal statistics literature (e.g., see Berliner,, 1996, for an early reference). This terminology is useful because it allows one to discuss highly parameterized joint distributions in terms of more manageable “pieces” (i.e., conditional and marginal distributions). The “data model” refers to the conditional distribution , where is an unobserved random variable assumed to be correlated across the indexes , , and , and will be used to denote a probability density function/probability mass function (pdf/pmf). Similarly, is referred to as the process model, where is a generic real-valued parameter vector. A parameter model, , is assumed for . Together, the “data model,” “process model,” and “parameter model” define the joint distribution that is used for statistical inference. That is,
[TABLE]
Each component of this hierarchical model has been discussed in the main text. Specifically, Section 2.1 describes the multinomial “data model,” the “process model” is defined in Sections 3.1 and 3.2, and the “parameter models are defined in Section 3.3. For a general review of the hierarchical modeling strategy see Cressie and Wikle, (2011).
Appendix A.i: Model Summary: The MN-STM
We now organize all the levels of the hierarchical model discussed in Sections 2.12.3 and 3.1 3.3. The MN-STM is proportional to the product of the following conditional and marginal distributions:
[TABLE]
where “ denotes a gamma density with shape and rate . To use (\thechapter), we need to specify , , , , , , , , , , , and . We choose these hyperparameters so that the prior is “flat.” Specifically, set (Gelman,, 2006). However, in general, these hyperparameters do not have to be the same, and one should consider alternate specifications. The full-conditional distributions for the shape parameters are computationally feasible to sample from using the adaptive rejection algorithm (Gilks and Wild,, 1992), which is possible due to Proposition 4. Although there are many levels in this statistical model, (\thechapter) is fairly parsimonious compared to the (latent Gaussian process model) MSTM. Both the MSTM and MN-STM includes parameters for regression parameters and random effects. The MN-STM only requires an additional shape parameters to estimate.
Using the collapsed Gibbs sampler we obtain replicates from the posterior distribution of , , and , which we denote with , , and , respectively; and . Then, the posterior replicate of is given by,
[TABLE]
Finally, to obtain posterior replicates of , first note that , so that . Then, we have
[TABLE]
which implies , and a posterior replicate can be found by computing . Continuing in this manner we obtain for , and . Finally, we compute averages and variances of (across ) to perform inference on .
One could also consider an empirical Bayesian version of (\thechapter) by substituting estimates of parameters and removing levels of the hierarchical model. Besides possibly simplifying computations, this approach may help address the issue of confounding without the need for the Moran’s I basis functions. In particular, empirical Bayesian hierarchical models treat as fixed, which may help with the issue of confounding.
Appendix A.ii: Basic Properties of the MN-STM Covariances
The parsimonious model presented in (\thechapter) allows for spatio-temporal dynamics, nonstationarity in space (see first paragraph of Section 3.3), nonstationarity in time, and asymmetric covariances. In this section, we clarify what is meant by spatio-temporal dynamics, non-stationary in time, and symmetric/asymmetric covariances.
We use the term “spatio-temporal dynamics” to refer to changes over time of the -dimensional vector,
[TABLE]
where and . Consider the “random walk” model as an example of a model without spatio-temporal dynamics. A random walk assumes with mutually independent of and . The conditional expected value is equal to an -dimensional vector of zeros.
The MN-STM incorporates spatio-temporal dynamics in the term , since for example,
[TABLE]
This expression is not identically equal to zero, and hence, we say that the MN-STM expresses spatio-temporal dynamics. Here the dynamics are determined by the propagator matrix .
Stationarity in time refers to the case where the covariance between multivariate spatial fields, and depend only on their temporal lag. We have that,
[TABLE]
which is not a function of , where we have assumed that . Stationarity is a special case, and occurs when and , so that
[TABLE]
for every and such that . However, our specifications in Section 2.3 imply that and for and , and hence, nonstationarity holds in time.
Symmetric covariances imply that the covariance between two spatial regions stay the same regardless of which variables/time points that are specified. In practice, the assumption of symmetry is unrealistic (see the discussion in Cressie and Wikle,, 2011, pg. 234). A popular class of symmetric covariances are known as “separable covariances” (e.g., see Stein,, 2005, and the references therein). Also, the covariances implied by a linear models for coregionalization can be interpreted as a type of mixture of separable covariancs, and also implies symmetry (e.g., see Jin et al.,, 2007, 2005, among others). The MN-STM yields asymmetric covariances. Specifically, the expression of the covariance (given ) between and is,
[TABLE]
The not equal to sign in the above holds because in our specification of the model in (\thechapter), , , , and .
Appendix A.iii: The MN-STM and the PMSTM
One has a choice between using the MN-STM or the PMSTM from Bradley et al., (2018b) in practice. To see this, change the inverse logit function with the following:
[TABLE]
Substituting this expression into the multinomial pdf we have,
[TABLE]
Then assuming we have,
[TABLE]
which is proportional to a Poisson distribution as a function of and the multivariate log-gamma distribution as a function of . Thus, one could use the PMSTM from Bradley et al., (2018b) to analyze in place of our MN-STM.
Appendix B: The Full-Conditional Distributions for the MN-STM
Each full conditional distribution associated with MN-STM is listed below.
The full conditional distribution for satisfies
[TABLE]
Rearranging terms we have
[TABLE]
which implies that is MLB with mean , covariance , shape , and , where
[TABLE]
and where in our implementation we set , , , and . Now, multiply the data model in (\thechapter) by , add to , and set the location parameter of in (\thechapter) to be equal to , where and is the orthonormal basis of the matrix immediately above. Let be the stacked vectors of each augmented vector (with improper prior) defined at the end of each step in this Gibbs sampler (not including ). Also let contain each shape parameter in (\thechapter). Then it follows from Proposition 2 that to simulate from , one can compute , where . 2. 2.
For and , the full conditional distribution for satisfies
[TABLE]
Rearranging terms we have
[TABLE]
which implies that is MLB with mean , covariance , shape , and , where
[TABLE]
and where in our implementation we set . Now, multiply the data model in (\thechapter) by , add to , add to the location parameter of in (\thechapter), and add to the location parameter of in (\thechapter), where and is the orthonormal basis of the matrix immediately above. Here, is , is , is , is , and is . Let be the stacked vectors of each augmented vector (with improper prior) defined at the end of each step in this Gibbs sampler (not including ). Then it follows from Proposition 2 that to simulate from , one can compute , where . 3. 3.
If , then the full conditional distribution for is given by
[TABLE]
Rearranging terms we have
[TABLE]
which implies that is MLB with mean , covariance , shape , and , where
[TABLE]
Now, multiply the data model in (\thechapter) by , add to , add to the location parameter of in (\thechapter), and add to the location parameter of in (\thechapter), where and is the orthonormal basis of the matrix immediately above. Here, is , is , is , is , and is . Let . Let be the stacked vectors of each augmented vector (with improper prior) defined at the end of each step in this Gibbs sampler (not including ). Then it follows from Proposition 2 that to simulate from , one can compute , where . 4. 4.
The full conditional distribution for satisfies
[TABLE]
Rearranging terms we have
[TABLE]
which implies that is MLB with mean , covariance , shape , and , where
[TABLE]
If then replace and with and within the expressions in (B.3). Now, multiply the data model in (\thechapter) by , add to , and add to the location parameter of in (\thechapter), where and is the orthonormal basis of the matrix immediately above. Let be the stacked vectors of each augmented vector (with improper prior) defined at the end of each step in this Gibbs sampler (not including ). Also let contain each shape parameter in (\thechapter). Then it follows from Proposition 2 that to simulate from , one can compute , where . 5. 5.
The full conditional distribution for is given by
[TABLE]
Rearranging terms we have
[TABLE]
which implies that is MLB with mean , covariance , shape , and , where
[TABLE]
Now, multiply the data model in (\thechapter) by , add to , and add to the location parameter of in (\thechapter), where and is the orthonormal basis of the matrix immediately above. Let be the stacked vectors of each augmented vector (with improper prior) defined at the end of each step in this Gibbs sampler (not including ). Also let contain each shape parameter in (\thechapter). Then it follows from Proposition 2 that to simulate from , one can compute , where . 6. 6.
The full conditional distribution for and , can be found as,
[TABLE]
where the -dimensional vector is defined in Steps 2 - 4 above. Take the first and second derivative of the log of (6) and obtain:
[TABLE]
It is well known that the log of the beta function is convex and the log of the gamma probability density function (with shape parameter greater than or equal to one) is concave (Dragomir et al.,, 2000). Thus, it follows from (B.8) that the pdf for is log-concave. The proofs for , , , , and are similar. 7. 7.
The full-conditional for was derived in Step 5, and can be found in (6). The derivation of the remaining shape parameters are straightforward and are as follows:
[TABLE]
[TABLE]
These full-conditional distributions are all computationally practical to simulate from using the adaptive rejection algorithm (see Step 5 for a proof of log-concavity).
Appendix C: A Latent Gaussian Process Model Representation of the Latent MLB model
Proposition 4 allows us to augment the model in (\thechapter) as follows:
[TABLE]
where is a Gaussian pdf with mean and covariance matrix , , , and . It follows from Proposition 3 that the latent MLB model in (\thechapter) is proportional to the density in (C.1) after integrating out the diagonal matrices , the diagonal matrices , and the diagonal matrix .
Acknowledgments
This research was partially supported by the U.S. National Science Foundation (NSF) and the U.S. Census Bureau under NSF grant SES-1132031, funded through the NSF-Census Research Network (NCRN) program. This article is released to inform interested parties of ongoing research and to encourage discussion. The views expressed on statistical issues are those of the authors and not necessarily those of the NSF or U.S. Census Bureau.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Abowd et al., (2013) Abowd, J., Schneider, M., and Vilhuber, L. (2013). “Differential privacy applications to Bayesian and linear mixed model estimation.” Journal of Privacy and Confidentiality , 5, 73–105.
- 2Abowd et al., (2009) Abowd, J., Stephens, B., Vilhuber, L., Andersson, F., Mc Kinney, K., Roemer, M., and Woodcock, S. (2009). “The LEHD infrastructure files and the creation of the Quarterly Workforce Indicators.” In Producer Dynamics: New Evidence from Micro Data , eds. T. Dunne, J. Jensen, and M. Roberts, 149–230. Chicago: University of Chicago Press for the National Bureau of Economic Research.
- 3Anderson, (1971) Anderson, T. (1971). The Statistical Analysis of Time Series . New York, NY: John Wiley and Sons.
- 4Andrieu et al., (2010) Andrieu, C., Doucet, A., and Holenstein, R. (2010). “Particle Markov chain Monte Carlo methods.” Journal of the Royal Statistical Society: Series B , 72, 269–342.
- 5Andrieu and Roberts, (2009) Andrieu, C. and Roberts, G. (2009). “The pseudo-marginal approach for efficient Monte Carlo computations.” The Annals of Statistics , 37, 697–725.
- 6Bakka et al., (2018) Bakka, H., Rue, H., Fuglstad, G.-A., Riebler, A., Bolin, D., Illian, J., Krainski, E., Simpson, D., and Lindgren, F. (2018). “Spatial modelling with R-INLA: A review.” ar Xiv preprint:1802.06350 .
- 7Banerjee et al., (2015) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2015). Hierarchical Modeling and Analysis for Spatial Data . London, UK: Chapman and Hall.
- 8Banerjee et al., (2008) Banerjee, S., Gelfand, A. E., Finley, A. O., and Sang, H. (2008). “Gaussian predictive process models for large spatial data sets.” Journal of the Royal Statistical Society, Series B , 70, 825–848.