Fitting spatial autoregressive logit and probit models using Stata: The spatbinary command

Abstract

Starting from version 15, Stata allows users to manage data and fit regressions accounting for spatial relationships through the sp commands. Spatial regressions can be estimated using the spregress, spxtregress, and spivregress commands. These commands allow users to fit spatial autoregressive models in cross-sectional and panel data. However, they are designed to estimate regressions with continuous dependent variables. Although binary spatial regressions are important in applied econometrics, they cannot be estimated in Stata. Therefore, I introduce spatbinary, a Stata command that allows users to fit spatial logit and probit models.

1 Introduction

In recent decades, spatial econometrics has been a growing field of study. In Stata software packages, spatial data were originally managed by a plethora of community contributed commands. Pisati (2001) and Drukker et al. (2013) released spatwmat and spmat, respectively, for managing spatial matrices. Other relevant contributions include software to perform spatial correlation tests (the spatcorr command of Pisati [2001]), geocode data (Ozimek and Miles 2011), calculate travel time (Huber and Rust 2016; Weber and Péclat 2017), and visualize detailed maps (Pisati 2018). Other than these utilities, community-contributed packages address spatial regression models in terms of cross-sectional data (Pisati 2001), panel data (Belotti, Hughes, and Mortari 2017), and endogenous regressors (Drukker, Prucha, and Raciborski 2013).

In addition, Stata 15 introduced sp, a suite of commands that allows users to manage spatial data and fit regressions with underlying spatial relationships. Some of the abovementioned community-contributed packages were embedded in sp.

As of Stata 17, official commands for spatial regression models (that is, spregress, spxtregress, and spivregress) are designed to fit models with continuous dependent variables. For models with binary responses, it is only possible to use spregress to fit spatial linear probability models. These models do not constrain the predicted values in [0, 1] (Wooldridge 2010). Conversely, spatial probit or logit models can be fit only by using other software. ¹

In applied economics, spatial models for binary outcomes have many uses. For instance, Klier and McMillen (2008) analyzed clustering in the automotive sector; Flores-Lagunes and Schnier (2012) used them to solve sample selection issues; Brasington, Flores-Lagunes, and Guci (2016) analyzed the adoption of open enrollment in schools; and more recently Prato et al. (2018) investigated the role of urban environment in pedestrian accidents.

The aim of this article is therefore to address the lack of Stata packages that are specifically aimed at fitting spatial models with binary response variables. Thus, I introduce spatbinary, a command that allows the estimation of logit and probit spatial autoregressive models. The spatbinary command is equivalent to spregress (with the option dvarlag) for binary dependent variables and to logit and probit for spatial data. The command is based on the generalized method of moments (gmm) estimators outlined in Pinkse and Slade (1998) and their approximations (Klier and McMillen 2008). The remainder of this article is organized as follows: the next section outlines the econometric setting for spatbinary, sections 3 and 4 present the syntax, section 5 presents some examples of the use of spatbinary, and section 6 concludes.

2 The econometric setting

Consider a binary outcome model in which the observed binary response y_i = 1 only if the unobserved latent variable u_i ∊ U is greater than 0; otherwise, y_i = 0. For instance, in the random utility setting (McFadden 1974), the latent variable represents the utility of a consumer facing the decision of whether to purchase a product or a service. If the underlying latent variable depends on a set of independent variables X and is also spatially autocorrelated, the binary spatial autoregressive (bsar) model assumes the form

$U=\rho WU+X\beta +\in$

1

where W is a row-standardized contiguity matrix, ∊ represents the error term, and parameters β and ρ are to be estimated. Spatial autocorrelation of the (latent) dependent variable U through parameter ρ implies a form of endogeneity because U appears in (1) both as a dependent variable and as a covariate (Arbia 2014). Specifically, this is a form of simultaneity that “arises when at least one of the explanatory variables is determined simultaneously along with the dependent variable” (Wooldridge 2010). In (1), it is straightforward to note that U is determined simultaneously with itself because it is the same variable.

Because U is unobserved, the spatial autocorrelation parameter ρ implies that the propensity of having a positive outcome is correlated with the propensity to have a positive outcome in nearby units of observation (Klier and McMillen 2008). A positive ρ implies clustering: units with a high probability of a positive outcome are located close to other peers characterized by a high probability that y_i = 1. Conversely, ρ < 0 means that the propensity is dispersed in space.

Depending on the distribution of the error term, (1) can be estimated through a logit or a probit model characterized by “transformed” independent variables X ^∗∗ , predicted probabilities P (y_i = 1), and generalized residuals e_i (table 1). The error-term variance is proportional to V = E( ∊ ^′ ∊ ) = [( I − ρ W ) ^′ ( I − ρ W )] ^{− 1} (Calabrese and Elkink 2014; McMillen 1992), and D contains the square root of the elements of V on its main diagonal.

Table 1.

Relevant quantities for logit and probit bsar models

	Probit	Logit
D	$\sqrt{\text{diag}(V)}$
X ∗∗	D − 1( I − ρ W ) − 1 X
P (yi = 1)	$\Phi \left(X_i^{**}\beta \right)$	$\frac{\exp \left(X_i^{**}\beta \right)}{1+\exp \left(X_i^{**}\beta \right)}$
ei	$\frac{\left\{y_i-P\left(y_i\right)\right\}\varphi \left(X_i^{**}\beta \right)}{P\left(y_i\right)*\left\{1-P\left(y_i\right)\right\}}$	yi − P (yi )

To estimate bsar parameters β and ρ, Pinkse and Slade (1998) proposed a gmm-type estimator (Hansen 1982) where coefficient estimates are chosen to minimize the quantity in (2).

$Q=N^{-}{\hspace{0.17em}}^1\left\{e{(\beta , \rho )}^{\prime }ZM{Z}^{\prime }e(\beta , \rho )\right\}$

2

Such a model requires a set of instruments Z that may include the independent variables X and their spatial lags WX ,…, W ⁿ X (Kelejian and Prucha 1998). Klier and McMillen (2008) noted that the estimator in (2) reduces to nonlinear two-stage least squares (n2sls) if M = ( Z ^′ Z ) ^{− 1} and proposed a simplified version of the n2sls model by introducing a linear approximation around a convenient starting point. The spatbinary command is designed to fit both models.

2.1 N2SLS estimator

The n2sls procedure for estimating model parameters Θ is the following:

1. Assume initial values Θ₀.

2. Repeat until convergence ² or until the maximum number of iterations is reached:

3. b. Obtain the fitted values ${\widehat{G}}_t$ from a linear regression of the gradient terms G _t = (∂ e _t )/(∂Θ _t ) on the instruments Z .

4. c. Update ${\Theta }_{t+1}={\Theta }_t+{\left[{\widehat{G}}_t^{\prime }{\widehat{G}}_t\right]}^{-1}{\widehat{G}}_t^{\prime }{e}_t$ .

The robust variance–covariance estimator is given by the following matrix:

${\left[{\widehat{G}}^{\prime }\widehat{G}\right]}^{-1}$ diag $\left(\widehat{G}{\text{'}}^{\prime }e{e}^{\prime }\widehat{G}\right){\left[\widehat{G}\text{'}\hspace{0.17em}\widehat{G}\right]}^{-1}$ .

2.2 Linearized estimator

The linearized estimator procedure is the following:

1. Through a nonspatial logit or probit model, obtain initial estimates of β ₀ and e ₀ (that is, assume ρ = 0).

2. Calculate the gradient G = (∂ e )/(∂Θ). These gradients are much simpler than in the n2sls case.

3. Regress G on Z and obtain the predicted gradient estimates $\widehat{G}$ .

4. Calculate $e_0+G_{\beta }^{\prime }{\beta }_0$ and regress it on $\widehat{G}$ . The obtained coefficients are the desired estimates of β and ρ.

2.3 Using the estimators in practice

The main advantage of the linearized model is computational. Indeed, it avoids the inversion of the matrix I −ρ W during each iteration. However, this advantage becomes less pronounced if the matrix W is small or sparse (that is, with many zero entries).

Furthermore, Klier and McMillen (2008) show that the linearized model provides a good approximation of the spatial autocorrelation parameter ρ only if the absolute value of the data-generating autocorrelation parameter lies in the interval 0.1–0.5. Outside this interval, estimates from the linearized model are upwardly biased (Arbia 2014). Moreover, assuming that the model is correctly specified, the standard errors produced by the linearized model are larger than those obtained by n2sls; thus, the linearized bsar is less efficient (Klier and McMillen 2008; Billé 2013). To summarize, the practical situations for which the linearized model can be used are the following:

• W is large.

• W is dense.

• The data have a weak spatial structure.

Otherwise, the n2sls approach is preferable, and coefficients from the linearized model can be used as starting values for the n2sls model. Neither model guarantees that the confidence intervals for ρ are bounded in the [−1, 1] interval (Billé 2013).

2.4 Measures of impact and marginal effects

As in standard probit and logit models, the coefficients β represent the impact of the independent variables X on the latent variable U. Such estimates are difficult to interpret because U is unobserved; however, they provide useful information about the direction of the impact of X on the probability that the outcome y is observed through their sign. To overcome this difficulty, a common strategy is to report measures of impact using margins (Williams 2012). In the general spatial autoregressive framework, a variation of an explanatory variable in a certain geographical unit affects the response variable both in the same location and in other locations because of spatial dependence (Anselin, Florax, and Rey 2004; Arbia 2014). Thus, measures of impact are split into direct and indirect components. The former measures a unit’s predicted contribution to its own probability of a positive outcome, and the latter measures the predicted impact of the other units’ contributions to a unit’s probability. In Stata, four types of measures of marginal impact are generally reported:

• the marginal effect (me), which represents the variation of the dependent variable in response to a unit variation in an explanatory variable (in margins, this is calculated using the dydx option);

• the elasticity, which represents the percent variation in the response variable in relation to a 1% variation of an explanatory variable (eyex option in margins);

• the semielasticity representing the percent variation in the response variable in relation to a unit variation of an explanatory variable (eydx option in margins); and

• a second type of semielasticity representing the variation in the dependent variable in natural units in relation to a 1% variation of an explanatory variable (dyex option in margins).

In the bsar framework, the response variable is P_i = P (y_i = 1). Table 2 summarizes the measures of impact related to the ith observation of covariate X ^k . In table 2, g_i is the first derivative of P_i with respect to the ith row of $X_i^k{\beta }_{X_i^k}$ . Such a derivative could be direct (that is, measuring a unit i’s predicted contribution to P_i ) or indirect (that is, measuring the other units’ contributions to the variation of P_i ). ³ In the spatbinary framework, the measures of impact in table 2 can be obtained using the postestimation command spatbinary_impact. Users can also use margins by exploiting the calculation of g_i using predict along with options directmargin, indirectmargin, and totalmargin.

Table 2.

Relevant measures of impact for bsar models

Measure	Expression
dydx	giβXk
eyex	$g_i{X}_i{\hspace{0.17em}}^k{\beta }_{X}k{[P_i]}^{-}{\hspace{0.17em}}^1$
eydx	giβXk[Pi ] − 1
dyex	$\begin{array}{l}{g}_i{X}_i{\hspace{0.17em}}^k{\beta }_{X}k\hfill \end{array}$

3 The spatbinary command

3.1 Syntax

Data should be spset before using spatbinary. The syntax of spatbinary is

spatbinary depvar indepvars [if] [in] [weight] , wmat( matname ) [ linearized n2sls logit probit noconstant force instr( varlist ) winstr( varlist ) impower( # ) noinstrconstant iterate( # ) start( matname ) tolerance( # ) display_options level( # ) coeflegend ]

fweight s are allowed; see [u] 11.1.6 weight.

3.2 Options

Main options

wmat( matname ) specifies that the spatial weight matrix W is stored in matrix matname. This option is required and is equivalent to dvarlag in spregress. The weight matrix must be created using spmatrix.

linearized requests that the linearized estimator be used. This is the default.

n2sls requests the n2sls estimator be used. Specifying both n2sls and linearized will cause spatbinary to ignore linearized.

logit requests to fit the logit bsar model. This is the default setting. If both probit and logit are specified, the program returns an error message and stops the execution of spatbinary.

probit requests to fit the probit bsar model.

noconstant suppresses the constant term in the model; see [r] Estimation options.

force requests that estimation be done when the estimation sample is a proper subset of the sample used to create the spatial weighting matrices. The default is to refuse to fit the model. Weighting matrices potentially connect all the spatial units. When the estimation sample is a subset of this space, the spatial connections differ, and spillover effects can be altered. In addition, the normalization of the weighting matrix differs from what it would have been had the matrix been normalized over the estimation sample. The better alternative to force is to first understand the spatial space of the estimation sample and then, if it is sensible, create new weighting matrices for it. See [sp] spmatrix and Missing values, dropped observations, and the W matrix in [sp] Intro 2.

Instruments options

The following options control the instruments specification. The number of instruments must be equal to or greater than the number of parameters to be estimated. Otherwise, spatbinary stops its execution.

instr( varlist ) specifies a list of instruments. By default, the independent variables are included.

winstr( varlist ) specifies a list of instruments to be spatially lagged using wmat( matname ). By default, the spatial lags of the independent variables are included. The option may be specified with instr( varlist ).

impower( # ) specifies the order of an instrumental-variable approximation used in fitting the model. The derivation of the estimator involves a product of # matrices. Increasing # may improve the precision of the estimation and will not cause harm but will require more computer time. The default is impower(1).

noinstrconstant omits the intercept from the instruments.

Options for the N2SLS model

These options affect the estimation only if n2sls is specified.

iterate( # ) specifies the number of iterations for n2sls estimation.

start( matname ) specifies that the starting values for the n2sls estimation be stored in vector matname. The default starting values are estimated using the linearized approach.

tolerance( # ) specifies the tolerance for the gmm criterion. The default is tolerance(1e-5).

Reporting

display_options: noci, nopvalues, noomitted, vsquish, noemptycells, baselevels, allbaselevels, nofvlabel, fvwrap( # ), fvwrapon( style ), cformat(% fmt ), pformat(% fmt ), sformat(% fmt ), and nolstretch; see [r] Estimation options.

level( # ) specifies the confidence level, as a percentage, for confidence intervals. The default is level(95) or as set by set level; see [u] 20.8 Specifying the width of confidence intervals.

coeflegend displays the legend instead of statistics.

3.3 Stored results

spatbinary stores the following in e()

4 Postestimation

After spatbinary, predict and spatbinary_impact are available. The latter is a wrapper of margins and estimates measures of impact such as marginal effects, elasticities, and semielasticities; it corresponds with estat impact following spregress; see [sp] spregress postestimation.

4.1 Syntax for predict

The available predictions are the following:

predict newvar , pr, the default, calculates the probability of a positive outcome.

predict newvar , xb calculates $X_i{\hspace{0.17em}}^{\ast \hspace{0.17em}\ast }\beta$ .

predict newvar , totalmargin calculates the total marginal effect with respect to $X_i{\hspace{0.17em}}^{\ast \hspace{0.17em}\ast }\beta$ . It is used with margins.

predict newvar , directmargin calculates the direct marginal effect of $X_i{\hspace{0.17em}}^{\ast \hspace{0.17em}\ast }\beta$ (g_i in section 2.4). It is used with margins.

predict newvar , indirectmargin calculates the indirect marginal effect of $X_i{\hspace{0.17em}}^{\ast \hspace{0.17em}\ast }\beta$ . It is used with margins.

predict newvar , residuals calculates probit or logit generalized residuals.

4.2 Syntax for spatbinary_impact

spatbinary_impact varlist [, [dydx| eyex| dyex| eydx] total direct indirect ]

Options dydx, eyex, dyex, and eydx are mutually exclusive. If total, direct, and indirect are all left unspecified, spatbinary_impact estimates them all.

4.2.1 Options

dydx calculates the marginal effect of varlist on the predicted probability (option pr in predict). dydx is the program default.

eyex calculates the elasticity of the predicted probability (option pr in predict) with respect to varlist.

dyex calculates the semielasticity (dP/dlnx) of the predicted probability (option pr in predict) with respect to varlist.

eydx calculates the semielasticity (dlnP/dx) of the predicted probability (option pr in predict) with respect to varlist.

total calculates the total measure of impact of varlist (defined by dydx, eyex, dyex, and eydx).

direct calculates the direct measure of impact. It captures own-unit contributions of varlist on a unit’s prediction.

indirect calculates the indirect measure of impact. It captures the contributions of the other units’ varlist on a unit’s prediction.

5 Examples

5.1 Setup

I illustrate spatbinary with five examples using homicide1990.dta (Messner et al. 2000); this dataset also provides examples for Stata’s official spatial regression command (spregress). Each observation in homicide1990 is a county in the southern United States.

Before we use spatbinary, the dataset must be declared to hold spatial data using spset, and a spatial weight matrix should be specified. To this end, homicide1990.dta comes with an ancillary shapefile storing the spatial information (coordinates and proximity) and is already spset. In the remainder of this section, the spatial weight matrix is created using spmatrix create and stored in W2, a row-standardized contiguity matrix. The examples include both probit and logit bsar models for illustrative purposes. Like their nonspatial counterparts, they generally give similar predictions (Greene 2018); however, the logit estimator is less likely to estimate probability estimates near 0 or 1. The choice between them depends on the specific application and field of research. Because coefficients are not readily interpretable, measures of impact are calculated. Thus, measures of impact throughout this section are related to gini.

The first three examples are devoted to logit and probit bsars with a dichotomized version of hrate, the county-level homicide rate per year per 100,000 persons, as the response variable. The independent variables are the Gini index, a measure of income inequality (gini), and the logarithm of the population (ln_population). Let us assume that the impact of income inequality on the probability of having a high homicide rate (also referred to as “the predicted probability” or “the probability that hrate_gt_p95=1”) is the main research interest of the analysis.

Examples 4 and 5 are based on simulated dependent variables Y1 and Y2. In both examples, the underlying latent variable is xb=-10+22*gini. Examples 4 and 5 are respectively based on ρ = 0.7 and ρ = 0.2.

5.2 Example 1—Linearized models with program defaults

This example illustrates the estimation of the linearized bsar using the default settings of spatbinary. In the first model, neither probit nor logit is specified; hence, a logit model is fit. The second model is a spatial probit. Further, because option n2sls is not specified, linearized probit and logit models are fit. Program defaults imply that the instruments are an intercept, the independent variables, and their spatial lag. This is equivalent to specifying options instr(ln_population gini), winstr(ln_population gini), and impower(1).

The coefficient rho is negative in both models and is significant only in the probit specification. Therefore, the propensity of having high homicide rates is characterized by weak dispersion: the probability of having hrate = 1 in a county decreases if the same propensity is high in the neighboring counties. The coefficient attached to gini is positive and highly significant in both models; this means that the higher the income inequality, the higher the probability of having high homicide rates. The coefficient of ln_population is positive and significant only in the probit model, suggesting a weak correlation between population and the probability of having high homicide rates. The package also provides Hansen’s test for overidentification, which applies when the number of instruments (moment conditions) is larger than the number of parameters to be estimated (see [r] gmm postestimation). In this example, there are four parameters to be estimated and five instruments. The test is not significant and safely suggests that the instruments are valid.

Because probit and logit are nonlinear models, their coefficients do not convey information about the magnitude of the relationship between the attached explanatory variable (gini and ln_population) and the probability that hrate_gt_p95=1. Moreover, the structure of the model [see (1) and section 2.4] implies spillovers through the existence of the spatial autocorrelation coefficient ρ. In this example, the implication is that the propensity of having high homicide rates is affected by own-county inequality (gini), own-county population (ln_population), and the propensity of having high homicide rates of neighboring counties, which in turn are affected by inequality and population. Thus, a county’s propensity of having high homicide rates is indirectly influenced by inequality and population in the neighboring counties. This fact adds complexity to the interpretation of the coefficients attached to the explanatory variables. mes are reported for the probit model because it exhibits a significant pattern of spatial autocorrelation (ρ is significant). mes for the logit model are also reported for appreciating the similarities between the two sets of me estimates.

In the output above, all the mes are significant. Furthermore, the total effect of gini as estimated by the probit model is equal to 1.31, meaning that, on average, an increase of 0.01 units of gini is related to an additional 1.31% probability of high homicide rates as a result of the summation of the direct and indirect effects. The direct me is equal to 1.84. ⁴ This finding suggests that an own-county 0.01-unit increase in the Gini index (gini) is to increase the probability of having high homicide rates by 1.84 percentage points. Furthermore, the indirect me is equal to −0.53. Hence, the across-county spillover effect of a 0.01-unit gini increase is associated with a reduction of probability of having high homicide rates by 0.53%. The negative sign of this me is coherent with the negative sign of ρ: the total effect of gini is positive, but the crosscounty component opposes to the own-county effect because of the negative spatial correlation (ρ = −0.375). Overall, the mes suggest that own-county income inequality increases the probability of having a high homicide rate, while other counties’ income inequality reduces the probability of having a high homicide rate.

As expected, the logit mes are remarkably similar to the corresponding probit me estimates. The total-effect point estimate is equal to 1.28 (1.31 in the probit model), and the direct effect is equal to 1.88 (1.84 in the probit model); lastly, the indirect me is equal to −0.60. However, the last is not significant; this reflects the fact that the nonsignificance of ρ in the logit model, and thus the increase of predicted probability, is likely associated only with own-county income inequality.

The same logit mes can be obtained by coding the expression shown in the first row of table 2 directly into the expression option of margins. This is recommended for the more advanced Stata users who want to exploit the flexibility of margins. The following Stata output uses margins to replicate the logit mes obtained above with spatbinary_impact. The reader is reminded that totalmargin, directmargin, and indirectmargin calculate the term g_i (see section 2.4).

5.3 Example 2—Beyond default instruments, the impower() option

In this example, the linearized models of example 1 are refit by replacing the default set of instruments with a more complex set. The task is achieved using the option impower(). A higher impower() value captures more complex spatial autocorrelation patterns and may result in more accurate estimates. For instance, Arbia (2014) suggests that impower(2) would eliminate the endogeneity implied in (1). Other contributions use impower(1)(Klier and McMillen 2008) or impower(3)(Calabrese and Elkink 2014). As in the previous example, leaving instr() and winstr() unspecified is equivalent to specifying instr(ln_population gini) and winstr(ln_population gini). In this example, impower() is equal to 3; the instruments are based on the independent variables that are multiplied by the first, second, and third powers of the spatial weight matrix W2. ⁵

Compared with example 1, the coefficients of gini and ln_population as well as the spatial correlation parameters are rather stable. Again, ρ is significant only in the probit model, and the test of overidentification restriction is insignificant. The total effect of gini on the predicted probability that hrate_gt_p95 = 1 is positive, and mes are reported to quantify such effects.

The gini me estimates for the probit model are very close to those of example 1. A 0.01-unit increase in gini is related to a 1.297% additional probability of high homicide rate. Again, the direct and indirect mes have an opposing effect. The former is equal to 1.792, and the latter is equal to −0.495. The low significance of the cross-county (indirect) me suggests that the dispersion in the propensity is weak; the effect of the explanatory variables is predominantly direct.

5.4 Example 3—N2SLS estimation

In this example, the n2sls model with default instruments is fit (that is, they are the same as example 1). Because the start() option is not specified, spatbinary internally fits the linearized model of example 1 to provide initial values. For both models, Hansen’s test for overidentification restrictions is safely insignificant, suggesting that the instruments are relevant.

The coefficient estimates are remarkably similar to those of the previous examples. However, rho is significant only at the 10% level of confidence, suggesting that the underlying spatial structure is weak. Moreover, the probit model standard errors of rho are larger than those reported in the corresponding linearized model (example 1): this may reflect model misspecification. Indeed, as reported in section 2.3, a correctly specified n2sls model should produce smaller standard errors than the corresponding linearized model.

Again, mes for the probit model are reported. Consistent with examples 1 and 2, the direct and total point estimates are positive, while indirect effects are negative. The total effect suggests that a 0.01-unit increase of income inequality as measured by the Gini index (gini) is associated with a 1.19% additional probability of having a high homicide rate. This is the combination of a positive direct effect and a negative crosscounty effect (−0.50, p = 0.05). In absolute value, mes are lower than those estimated in examples 1 and 2 by the linearized models.

5.5 Example 4—Simulated data with ρ = 0.7

In this example, the models are correctly specified because they are based on the datagenerating function. The first bsar logit model is based on the linearized estimator, and the second and third are based on n2sls. The estimated value of ρ is maximum in the first regression. This is expected because the data-generating value is higher than 0.5 (the linearized model is upwardly biased). However, the real value (0.7) is included in the confidence interval. In general, the true coefficients (−10 for the intercept and 22 for gini) are safely included in their related confidence intervals.

The n2sls models have lower standard errors; this was expected because this estimator is more efficient than the linearized version. In the second model, the number of instruments (3) is equal to the number of parameters to be estimated; hence, instrument validity cannot be tested. To perform this test, impower(2) is specified in the third regression. Specifically, Hansen’s test in the third model implies that the second-order instruments (a constant, gini, W2*gini, and (W2^2)*gini) are valid. Overall, the positive sign of ρ suggests that the propensities of Y1=1 are clustered in space: counties with high probability of Y1=1 are located close to each other. This pattern also suggests that the cross-county (indirect) effect of the explanatory variables have the same sign of the direct effect.

To give more precise quantitative indications about the effect of gini on the predicted probability that Y1=1, I report the four types of measures of impact (see section 2.4). As the choice between the logit and probit models, the type of measure of impact to be reported depends on the context of research. Measures of impact are estimated using margins for illustrative purposes. In the following Stata outputs, the order of reporting is total impact, direct impact, and indirect impact. The following command calculates the mes of gini. It is equivalent to spatbinary_impact gini, dydx. The mes depict that a 0.01-unit increase in the Gini index yields to a 2.31% additional predicted probability that Y1=1. Both the direct and indirect effects are positive; the larger component is related to the cross-county effect of gini because the indirect effect is equal to 1.43 and the direct effect is equal to 0.88.

Elasticity estimates are reported below. Specifying spatbinary_impact gini, eyex gives the same results. The elasticity represents the proportional increase of the predicted probability associated with a 1% proportional (that is, multiplicative) increase of gini. The total elasticity is equal to 24. Again, the indirect effect (15.11) is stronger than the direct effect (9.38).

The semielasticities giving the additive increase in probability that Y1=1 in response to a proportional increase of gini are reported in the Stata output below. They can be obtained by specifying spatbinary_impact gini, dyex. This means that a 1% (multiplicative) increase of gini leads to approximately a percentage point of the predicted probability that Y1=1. Specifically, this increase is mostly attributable to the cross-county (indirect) effect: a 1% proportional increase of gini is related to a 0.62% positive variation in terms of predicted probability.

Lastly, to calculate the semielasticities giving the proportional increase in probability of Y1=1 in response to a unit increase of gini, Stata users can plug into margins the expression predict(totalmargin)*_b[Y1:gini]/predict(pr). The equivalent postestimation command is spatbinary_impact gini, eydx. As expected, the indirect effect is stronger than the direct effect. In quantitative terms, such an effect implies that a 0.01-unit increase in gini in other geographical units is related to a 38.7% proportional increase in the predicted probability.

5.6 Example 5—Simulated data with ρ = 0.2

As in the previous example, the models are correctly specified. Because ρ is small, the linearized model estimates and the first set of n2sls coefficients are essentially the same, with a gain in efficiency related to the standard errors of the independent variables. Again, the data-generating coefficients are safely included in their related confidence intervals. The models highlight that the total effect of gini on the propensity of Y1=1 is positive. Clustering of such propensity (ρ > 0) implies that the direct and indirect effects of gini act in the same direction.

6 Conclusions

This article discussed the implementation of spatbinary, a Stata package that allows users to fit bsar logit and probit models. The spatbinary command extends the suite of community-contributed and official commands for spatial regression. A possible development of spatbinary would be its extension to multinomial outcomes.

Supplemental Material

Supplemental Material, sj-zip-1-stj-10.1177_1536867X221106373 - Fitting spatial autoregressive logit and probit models using Stata: The spatbinary command

Supplemental Material, sj-zip-1-stj-10.1177_1536867X221106373 for Fitting spatial autoregressive logit and probit models using Stata: The spatbinary command by Daniele Spinelli in The Stata Journal