CONTRIBUTED RESEARCH ARTICLES 130 Estimating Spatial Probit Models in R

In this article we present the Bayesian estimation of spatial probit models in R and provide an implementation in the package spatialprobit. We show that large probit models can be estimated with sparse matrix representations and Gibbs sampling of a truncated multivariate normal distribution with the precision matrix. We present three examples and point to ways to achieve further performance gains through parallelization of the Markov Chain Monte Carlo approach.


Introduction
The abundance of geolocation data and social network data has lead to a growing interest in spatial econometric methods, which model a contemporaneous dependence structure using neighboring observations (e.g.friends in social networks).
There are a variety of R packages for spatial models available, an incomplete list includes spBayes (Finley and Banerjee, 2013), spatial (Venables and Ripley, 2002), geostatistical packages called geoR (Diggle and Ribeiro, 2007) and sgeostat (Majure and Gebhardt, 2013), spdep (Bivand, 2013), sphet (Piras, 2010), sna (Butts, 2013) and network (Butts et al., 2013).While all of these packages deal with linear spatial models, in this article we focus on a nonlinear model, the spatial probit model, and present the Bayesian estimation first proposed by LeSage (2000).As will be shown below, one crucial point we have been working on was the generation of random numbers of a truncated multivariate normal distribution in very high dimensions.
Together with these high dimensions and the size of problems in spatial models comes the need to work with sparse matrix representations rather than the usual dense matrix().Popular R packages for dealing with sparse matrices are Matrix (Bates and Maechler, 2013) and sparseM (Koenker and Ng, 2013).
In the next section we first introduce spatial probit models and describe the Bayesian estimation procedure of the SAR probit model in detail.Subsequently we discuss some implementation issues in R, describe the sarprobit method in package spatialprobit (Wilhelm and de Matos, 2013), compare it against maximum likelihood (ML) and generalized method of moments (GMM) estimation in McSpatial (McMillen, 2013) and illustrate the estimation with an example from social networks.We also illustrate how to parallelize the Bayesian estimation with the parallel package.

Spatial probit models
The book of LeSage and Pace (2009) is a good starting point and reference for spatial econometric models in general and for limited dependent variable spatial models in particular (chapter 10, p. 279).
Suppose we have the spatial autoregressive model (SAR model, spatial lag model) for z = (z 1 , . . ., z n ) with some fixed matrix of covariates X (n × k) associated with the parameter vector β (k × 1).The matrix W (n × n) is called the spatial weight matrix and captures the dependence structure between neighboring observations such as friends or nearby locations.The term Wz is a linear combination of neighboring observations.The scalar ρ is the dependence parameter and will assumed abs(ρ) < 1.The k + 1 model parameters to be estimated are the parameter vector β and the scalar ρ.
In a spatial probit model, z is regarded as a latent variable, which cannot be observed.Instead, the observables are only binary variables y i (0, 1) as y i can reflect any binary outcome such as survival, a buy/don't buy decision or a class variable in binary classification problems.For identification, σ 2 is often set to σ 2 = 1.
The R Journal Vol.5/1, June ISSN 2073-4859 The data generating process for z is Note that if ρ = 0 or W = I n , the model reduces to an ordinary probit model, for which Wooldridge (2002, chapter 17) and Albert (2007, section 10.1) are good references.The ordinary probit model can be estimated in R by glm().
Another popular spatial model is the spatial error model (SEM) which takes the form where as before, only binary variables y i (0,1) can be observed instead of z i .Our following discussion focuses on the spatial lag probit model ( 1), but the estimation of the probit model with spatial errors is covered in spatialprobit as well.

Bayesian estimation
Maximum Likelihood and GMM estimation of the spatial probit is implemented in the package McSpatial (McMillen, 2013) with the methods spprobitml and spprobit.Here we present the details of Bayesian estimation.Although GMM is massively more efficient computationally than Bayesian Markov Chain Monte Carlo (MCMC), as a instrumental-variables estimation it will work well only in very large samples.This point is also brought by Franzese et al. (2013), who compares different spatial probit estimation strategies.
The basic idea in Bayesian estimation is to sample from a posterior distribution of the model parameters p(z, β, ρ|y) given the data y and some prior distributions p(z), p(β), p(ρ).See for example Albert (2007) and the accompanying package LearnBayes for an introduction to Bayesian statistics in R (Albert, 2012).This sampling for the posterior distribution p(z, β, ρ|y) can be realized by a Markov Chain Monte Carlo and Gibbs sampling scheme, where we sample from the following three conditional densities p(z|β, ρ, y), p(β|z, ρ, y) and p(ρ|z, β, y): 1. Given the observed variables y and parameters β and ρ, we have p(z|β, ρ, y) as a truncated multinormal distribution subject to z i ≥ 0 for y i = 1 and z i < 0 for y i = 0, which can be efficiently sampled from using the method rtmvnorm.sparseMatrix() in package tmvtnorm (Wilhelm and Manjunath, 2013), see the next section for details.
2. For a normal prior β ∼ N(c, T), we can sample p(β|ρ, z, y) from a multivariate normal as The standard way for sampling from this distribution is rmvnorm() from package mvtnorm (Genz et al., 2013).
3. The remaining conditional density p(ρ|β, z, y) is which can be sampled from using Metropolis-Hastings or some other sampling scheme (i.e.Importance Sampling).We implement a grid-based evaluation and numerical integration proposed by LeSage and Pace and then draw from the inverse distribution function (LeSage and Pace, 2009, p. 132).

Implementation issues
Since MCMC methods are widely considered to be (very) slow, we devote this section to the discussion of some implementation issues in R. Key factors to estimate large spatial probit models in R include the usage of sparse matrices and compiled Fortran code, and possibly also parallelization, which has been introduced to R 2.14.0 with the package parallel.We report estimation times and memory requirements for several medium-size and large-size problems and compare our times to those reported by LeSage and Pace (2009, p. 291).We show that our approach allows the estimation of large spatial probit models in R within reasonable time.
Drawing p(z|β, ρ, y) In this section we describe the efforts made to generate samples from p(z|β, ρ, y) subject to z i ≥ 0 for y i = 1 and z i < 0 for y i = 0.The generation of samples from a truncated multivariate normal distribution in high dimensions is typically done with a Gibbs sampler rather than with other techniques such as rejection sampling (for a comparison, see Wilhelm and Manjunath (2010)).The Gibbs sampler draws from univariate conditional densities f (z The distribution of z i conditional on z −i is univariate truncated normal with variance and mean and bounds −∞ ≤ z i ≤ 0 for y i = 0 and 0 ≤ z i ≤ ∞ for y i = 1. Package tmvtnorm has such a Gibbs sampler in place since version 0.9, based on the covariance matrix Σ.Following Geweke (1991Geweke ( , 2005, p. 171), p. 171), the Gibbs sampler is easier to state in terms of the precision matrix H, simply because it requires fewer and easier operations in the above equations ( 7) and ( 9).These equations will further simplify in the case of a sparse precision matrix H, in which most of the elements are zero.With sparse H, operations involving H i,−i need only to be executed for all non-zero elements rather than for all n − 1 variables in z −i .
In our spatial probit model, the covariance matrix Σ = [(I n − ρW) (I n − ρW)] −1 is a dense matrix, whereas the corresponding precision matrix H = Σ −1 = (I n − ρW) (I n − ρW) is sparse.Hence, using Σ for sampling z is inefficient.For this reason we reimplemented the Gibbs sampler with the precision matrix H in package tmvtnorm instead (Wilhelm and Manjunath, 2013).
Suppose one wants to draw N truncated multinormal samples in n dimensions, where the precision matrix H is sparse (n × n) with only m < n entries different from zero per row on average (e.g.m = 6 nearest neighbors or average branching factor in a network).With growing dimension n, two types of problems arise with a usual dense matrix representation for H: 1. Data storage problem: Since matrix H is n × n, the space required to store the dense matrix will be quadratic in n.One remedy is, of course, using some sparse matrix representation for H (e.g.packages Matrix, sparseM etc.), which actually holds only n • m elements instead of n • n.
The following code example shows the difference in object sizes for a dense vs. sparse identity matrix I n for n = 10000 and m = 1.

Mb
The R Journal Vol.rtmvnorm.sparseMatrix()(tmvtnorm) for a varying number of samples N and dimension n.The precision matrix in each case is H = (I n − ρW) (I n − ρW), the spatial weight matrix W contains m non-zero entries per row.Times are in seconds and measured on an Intel® Core™i7-2600 CPU @3.40 GHz.
One more performance issue we discuss here is the burn-in size in the innermost Gibbs sampler generating z from p(z|β, ρ, y).Depending on the start value, Gibbs samplers often require a certain burn-in phase until the sampler mixes well and draws from the desired target distribution.In our MCMC setup, we only draw one sample of z from p(z|β, ρ, y) in each MCMC iteration, but possibly in very high dimensions (e.g.N = 1, n = 10000).With burn.in=2samples, we have to generate 21 draws in order to keep just one.In our situation, a large burn-in size will dramatically degrade the MCMC performance, so the number of burn-in samples for generating z has to be chosen carefully.LeSage and Pace (2009) discuss the role of the burn-in size and often use burn.in=1, when the Gibbs sampler starts from zero.Alternatively, they also propose to use no burn-in phase at all (e.g.burn.in= ), but then to set the start value of the sampler to the previous value of z instead of zero.

QR decomposition of (I n − ρW)
The mean vector µ of the truncated normal samples z in equation (3) takes the form However, inverting the sparse matrix S = (I n − ρW) will produce a dense matrix and will therefore offset all benefits from using sparse matrices (i.e.memory consumption and size of problems that can be solved).It is preferable to determine µ by solving the equations (I n − ρW)µ = Xβ with a QR decomposition of S = (I n − ρW).We point out that there is a significant performance difference between the usual code mu <-qr.solve(S,X%*% beta) and mu <-solve(qr(S),X %*% beta).The latter will apply a QR decomposition for a sparse matrix S and will use a method qr() from the package Matrix, whereas the first function from the base package will coerce S into a dense matrix.

Drawing p(β|z, ρ, y)
Drawing p(β|z, ρ, y) in equation ( 4) is a multivariate normal distribution, whose variance T * does not depend on z and ρ.Therefore, we can efficiently vectorize and generate temporary draws β tmp ∼ N(0, T * ) for all MCMC iterations before running the chain and then just shift these temporary draws by c * in every iteration to obtain a β ∼ N(c * , T * ).
Determining log-determinants and drawing p(ρ|z, β, y) The computation of log-determinants for ln |I n − ρW|, whose evaluation is frequently needed for drawing p(ρ|z, β, y), becomes a challenging task for large matrices.Bivand (2010) gives a survey of the different methods for calculating the Jacobian like the Pace and Barry (1997) grid evaluation of ρ in the interval [−1, . . ., +1] and spline approximations between grid points or the Chebyshev approximation of the log-determinant (Pace and LeSage, 2004).All of them are implemented in package spdep.For the estimation of the probit models we are using the existing facilities in spdep (do_ldet()).

Computation of marginal effects
In spatial lag models (and its probit variants), a change of an explanatory variable x ir will affect both the same response variable z i (direct effects) and possibly all other responses z j (j = i; indirect effects or spatial spillovers).The resulting effects matrix S r (W) is n × n for x r (r = 1, . . ., k) and three summary measures of S r (W) will be computed: average direct effects (M r (D) = n −1 trS r (W)), average total effects (M r (T) = n −1 1 n S r (W)1 n ) and the average indirect effects as the difference of total and direct effects (M r (I) = M r (T) − M r (D)).In contrast to the SAR probit, there are no spatial spill-over effects in the SEM probit model (2).In the MATLAB spatial econometrics toolbox (LeSage, 2010), the computation of the average total effects requires an inversion of the n × n matrix S = (I n − ρW) at every MCMC iteration.Using the same QR-decomposition of S as described above and solving the equation solve(qr(S),rep(1,n)) speeds up the calculation of total effects by magnitudes.

Examples Package spatialprobit
After describing the implementation issues to ensure that the estimation is fast enough and capable to handle large problems, we briefly describe the interface of the methods in package spatialprobit (Wilhelm and de Matos, 2013) and then turn to some examples.
The main estimation method for the SAR probit model sar_probit_mcmc(y,X,W) takes a vector of dependent variables y, a model matrix X and a spatial weight matrix W. The method sarprobit(formula, W,data) is a wrapper which allows a model formula and a data frame.Both methods require a spatial weight matrix W to be passed as an argument.Additionally, the number of MCMC start values, the number of burn-in iterations and a thinning parameter can be specified.The estimation fit <-sarprobit(y ~x,W,data) returns an object of class sarprobit.The model coefficients can be extracted via coef(fit).summary(fit) returns a coefficient table with z-values, impacts(fit) gives the marginal effects and the plot(fit) method provides MCMC trace plots, posterior density plots as well as autocorrelation plots of the model parameters.logLik(fit) and AIC(fit) return the log likelihood and the AIC of the model for model comparison and testing.

Experiment from LeSage/Pace
We replicate the experiment from LeSage and Pace (2009, section 10.1.5)for n = 400 and n = 1000 random points in a plane and a spatial weight matrix with the 6 nearest neighbors.The spatial probit model parameters are σ 2 = 1, β = (0, 1, −1) and ρ = 0.75.We generate data from this model with the following code.
Table 3 shows that a change of the Gibbs sampler burn-in size m for drawing z has little effect on the mean and the variance of the estimates.m = 1 means taking the previous value of z, m = 2, 5, 10 start the chain from z = 0. Clearly, choosing m = 2 might not be enough as burn-in phase.

Random graph example
As a second example we present the estimation of the probit model based on a random undirected graph with n = 200 nodes and an average branching factor of 3. Of course this estimation procedure can be applied to real network structures such as from social networks.The package igraph can be used to create random graphs and to compute the adjacency matrix A as well as the spatial weight matrix W (Csardi and Nepusz, 2006).

Diffusion of innovation and information
A last example looks at the information flow in social networks.Coleman et al. (1957Coleman et al. ( , 1966) have studied how innovations (i.e.new drugs) diffuse among physicians, until the new drug is widely adopted by all physicians.They interviewed 246 physicians in 4 different US cities and looked at the role of 3 interpersonal networks (friends, colleagues and discussion circles) in the diffusion process (November, 1953-February, 1955).Figure 4 illustrates one of the 3 social structures based on the question "To whom do you most often turn for advice and information?".The data set is available at http://moreno.ss.uci.edu/data.html#ckm,but also part of spatialprobit (data(CKM)).See also Burt (1987) and den Bulte and Lilien (2001) for further discussion of this data set.The dependent variable in the model is the month in which a doctor first prescribed the new drug.Explanatory variables are the social structure and individual variables.
The R Journal Vol.5/1, June ISSN 2073-4859 The R Journal Vol.We are estimating a standard probit model as well as a SAR probit model based on the "advisorship" network and trying to find determinants for all adopters of the new drug.We do not aim to fully reanalyze the Coleman data set, nor to provide a detailed discussion of the results.We rather want to illustrate the model estimation in R. We find a positive relationship for the social influence variable in the probit model, but social contagion effects as captured by ρW in the more sound SAR probit model is rather small and insignificant.This result suggests that social influence is a factor in information diffusion, but the information flow might not be correctly described by a SAR model.Other drivers for adoption which are ignored here, such as marketing efforts or aggressive pricing of the new drug may play a role in the diffusion process too (den Bulte and Lilien, 2001).

Parallel estimation of models
MCMC is, similar to the bootstrap, an embarrassingly parallel problem.It can be easily run in parallel on several cores.From version 2.14.0,R offers a unified way of doing parallelization with the parallel package.There are several different approaches available to achieve parallelization and not all approaches are available for all platforms.See for example the conceptual differences between the two main methods mclapply and parLapply, where the first will only work serially on Windows.Users are therefore encouraged to read the parallel package documentation for choosing the appropriate way.

Summary
In this article we presented the estimation of spatial probit models in R and pointed to the critical implementation issues.Our performance studies showed that even large problems with n = 10, 000 or n = 100, 000 observations can be handled within reasonable time.We provided an update of tmvtnorm and a new package spatialprobit on CRAN with the methods for estimating spatial probit models implemented (Wilhelm and de Matos, 2013).The package currently implements three limited dependent models: the spatial lag probit model (sarprobit()), the probit model with spatial errors (semprobit()) and the SAR Tobit model (sartobit()).The Bayesian approach can be further extended to other limited dependent spatial models, such as ordered probit or models with multiple spatial weights matrices.We are planning to include these in the package in near future.

Figure 1 :
Figure 1: Random graph with n = 200 nodes and average branching factor 3

Figure 2 :Figure 3 :Figure 4 :
Figure 2: MCMC trace plots for model parameters with horizontal lines marking the true parameters

Table 1 :
Even with a sparse matrix representation for H, a naive strategy that tries to access an arbitrary matrix element H[i,j] like in triplet or hash table representations, will result in N • n • n matrix accesses to H.This number grows quadratically in n.For example, for N = 30 and n = 10000 it adds up to 30 • 10000 • 10000 = 3 billion hash function calls, which is inefficient and too slow for most applications.Iterating only over the non-zero elements in row i, for example in term H i,−i , reduces the number of accesses to H to only N • n • m instead of N • n • n, which furthermore will only grow linearly in n for a fixed m.Suitable data structures to access all non-zero elements of the precision matrix H are linked lists of elements for each row i (list-of-lists) or, thanks to the symmetry of H, the compressed-sparse-column/row representation of sparse matrices, directly available from the Matrix package.How fast is the random number generation for z now?We performed a series of tests for varying sizes N and n and two different values of m.The results presented in Table1show that the sampler with sparse H is indeed very fast and works fine in high dimensions.The time required scales with N • n.Performance test results for the generation of truncated multivariate normal samples z with 5/1, June ISSN 2073-4859 2. Data access problem: - ***' . 1 '**' . 1 '*' .5 '.' .1 ' ' 1 Figures 2 and 3 show the trace plots and posterior densities as part of the MCMC estimation results.Table4compares the SAR probit and standard probit estimates.

Table 2 :
Bayesian SAR probit estimates for n = 400 and n = 1000 with N = 1000 draws and 200 burn-in samples.m is the burn-in size in the Gibbs sampler drawing z.Timings in R include the computation of marginal effects and were measured using R 2.15.2 on an Intel® Core™i7-2600 CPU @3.40 GHz.Times marked with ( †) are taken from LeSage and Pace (2009, Table10.1).To allow for a better comparison, these models were estimated anew ( ‡) using MATLAB R2007b and the spatial econometrics toolbox (LeSage, 2010) on the very same machine used for getting the R timings.The ML and GMM timings ( §) do not involve the computation of marginal effects.

Table 3 :
Effects of the Gibbs sampler burn-in size m on SAR probit estimates for n = 400, N = 1000 draws and burn.in=2

Table 4 :
SAR probit estimates vs. probit estimates for the random graph example

Table 5 :
SAR probit and standard probit estimates for the Coleman data.