# Land-Use Filtering for Nonstationary Spatial Prediction of Collective
Efficacy in an Urban Environment

J. Brandon<EMAIL_ADDRESS>[ Christopher
R<EMAIL_ADDRESS>[ Bethany
<EMAIL_ADDRESS>[ Nicolo
<EMAIL_ADDRESS>[ Catherine A.
<EMAIL_ADDRESS>[ Department of
Statistics and Data Sciences, University of Texas at Austinpresep=, ]e1,e2
Department of Sociology, The Ohio State Universitypresep=, ]e4,e3  Population
Research Institute, The Ohio State Universitypresep=, ]e5

###### Abstract

Collective efficacy – the capacity of communities to exert social control
toward the realization of their shared goals – is a foundational concept in
the urban sociology and neighborhood effects literature. Traditionally,
empirical studies of collective efficacy use large sample surveys to estimate
collective efficacy of different neighborhoods within an urban setting. Such
studies have demonstrated an association between collective efficacy and local
variation in community violence, educational achievement, and health. Unlike
traditional collective efficacy measurement strategies, the Adolescent Health
and Development in Context (AHDC) Study implemented a new approach, obtaining
spatially-referenced, place-based ratings of collective efficacy from a
representative sample of individuals residing in Columbus, OH. In this paper,
we introduce a novel nonstationary spatial model for interpolation of the AHDC
collective efficacy ratings across the study area which leverages
administrative data on land use. Our constructive model specification strategy
involves dimension expansion of a latent spatial process and the use of a
filter defined by the land-use partition of the study region to connect the
latent multivariate spatial process to the observed ordinal ratings of
collective efficacy. Careful consideration is given to the issues of parameter
identifiability, computational efficiency of an MCMC algorithm for model
fitting, and fine-scale spatial prediction of collective efficacy.

Bayesian statistics,

Data augmentation,

Dimension expansion,

Nonstationarity,

Sociology,

Spatial statistics,

###### keywords:

, , , , and

##  1 Introduction

For decades, social science research has investigated how neighborhood
residents can mobilize to address local problems (Shaw and McKay, 1942;
Jacobs, 1961). The theory of neighborhood collective efficacy has been highly
influential in this respect (Sampson, Raudenbush and Earls, 1997). Uniting
Coleman’s 1988 concept of social capital and Bandura’s 1982; 1986 research on
personal and group-level self-efficacy, neighborhood collective efficacy is a
construct capturing the collective capacity of community members to exert
social control toward the realization of their shared goals (e.g., low levels
of crime).

Empirical tests of the neighborhood collective efficacy theory often rely on
data collected through sample surveys administered to residents of different
neighborhoods within an urban center. Collective efficacy survey instruments
vary, but typically ask urban residents to report perceptions of local trust,
monitoring, and intervention norms in their neighborhood. We refer to these
three categories of questions – ”trust,” ”observation,” and ”defense,”
respectively – as the components of collective efficacy. Reports on these
components are then aggregated within geographic units consistent with the
notion of neighborhoods (e.g., census tracts). For example, the Project on
Human Development in Chicago Neighborhoods (PHDCN), conducted in the mid
1990s, used neighborhood ratings from nearly 8,000 Chicago residents clustered
within 343 census-based neighborhoods to capture collective efficacy
Raudenbush and Sampson, 1999; see also Matsueda and Drakulich, 2016; Wickes et
al., 2019. In the PHDCN and other studies, multilevel/hierarchical regression
models were used to capture between-neighborhood variation in collective
efficacy through random effects, which are assumed to be spatially independent
or dependent over the study region depending on the study design. Studies of
neighborhood-level collective efficacy have found inverse associations with a
host of neighborhood problems such as rates of homicide (Sampson, Raudenbush
and Earls, 1997), intimate partner violence (Browning, 2002), child
maltreatment (Molnar et al., 2015), and chronic disease (Cohen et al., 2006),
even after adjusting for neighborhood sociodemographic factors.

In this paper, we extend the notion of collective efficacy as a continuously-
indexed spatial feature of an urban area, which is defined by the collective
impressions of places individuals visit as part of their normal, everyday
routine. Using point-referenced collective efficacy data collected as part of
the Adolescent Health and Development in Context (AHDC) Study, described in
detail in Section 2, we estimate the components of collective efficacy across
the city of Columbus, OH. Our approach also uses taxation-based records on the
land use of parcels – small geographic units often used as a proxy for a
“place.” Consistent with recent work suggesting that land-use compositions
affect the levels of collective efficacy within neighborhoods (Corcoran et
al., 2018), we use the category of land use to model spatial variation in the
mean of the collective efficacy component processes. Our novel contribution
from a statistical perspective is that we capture land-use driven (i.e. driven
by spatially-referenced covariate information), second-order nonstationarity
through a dimension-expansion strategy, building upon existing approaches for
covariate-driven nonstationary spatial modeling (Calder, 2008; Schmidt,
Guttorp and O’Hagan, 2011; Reich et al., 2011; Neto, Schmidt and Guttorp,
2014; Ingebrigtsen, Lindgren and Steinsland, 2014; Risser and Calder, 2015;
Risser et al., 2019). In a similar spirit to the dimension expansion approach
to nonstationary spatial modeling of Bornn, Shaddick and Zidek (2012), we
augment the dimension of the latent collective efficacy component processes to
the number of land-use categories and model a multivariate, stationary spatial
process as if spatially-referenced data on land-use specific collective
efficacy components were observable everywhere. We then relate the
observations to the higher dimensional random process using a filter defined
by the land-use partitioned study area. We show that our land-use filtering
model provides a better fit to the AHDC collective efficacy data than a
traditional spatial generalized linear mixed model with a single latent
second-order stationary Gaussian process (Banerjee, Carlin and Gelfand, 2014).

Our proposed methodology advances the measurement of collective efficacy in
two ways. First, it is compatible with the more cost effective data collection
strategy employed in the AHDC Study, where study participants report on
collective efficacy levels at their routine activity locations, of which most
have between five and eight. In addition, our model readily allows point-level
prediction of the components of collective efficacy across the entire study
area to better understand within neighborhood variation in collective
efficacy. Point-level variability in collective efficacy may help explain
phenomena such as crime which has been shown to concentrate at particular
locations within neighborhoods (Weisburd et al., 2016). Finally, we note that
while developed for the study of collective efficacy, our methodology can be
readily applied in other spatial prediction settings where the study area can
be partitioned into an arbitrary number of land-use categories.

The outline of the paper is as follows. In Section 2, we introduce the AHDC
Study and the collective efficacy data. In addition, we describe exploratory
analyses of the data that motivate our novel land-use filtering methodology.
Section 3 introduces our land-use filtering model and details land-use
filtering as applied to our ordinal response variable. We provide a summary of
inferences and predictions from our fitted model to the AHDC data in Section 4
and explore model performance when data are generated from a land-use
filtering process through a simulation study in Section 5. Lastly, we discuss
implications of our modeling and data collection choices in Section 6.

##  2 Data and exploratory analyses

In this section, we describe the AHDC data in more detail and visually
summarize the spatial patterning of the ratings of the components of
collective efficacy. We also introduce the land-use data, which drive our
proposed spatial filtering strategy for smoothing the ratings across the study
area. This section concludes with summaries of preliminary models that
motivate the more complex modeling strategy introduced in Section 3.

###  2.1 AHDC data

The AHDC Study is a longitudinal data collection project designed to improve
understanding of how social, psychological, and biological processes shape
youth developmental outcomes. Participants in the study are members of a
representative sample of households with youth aged 11 to 17 residing within
the I-270 belt loop in Franklin County, OH, which contains the city of
Columbus and some of its interior suburbs. In this paper, we use data from
AHDC Wave 1, which was collected between 2014-2016. For each sampled
household, one randomly-sampled youth aged 11-17 was selected for
participation in the study. Upon enrollment, the primary caregiver, typically,
but not always, the mother of the enrolled youth, filled out an entrance
questionnaire which included questions on family and household composition,
alcohol and substance use, employment and income, health, and social support.
In addition to answering these demographic and socioeconomic background
questions, caregivers also listed locations they frequent as part of their
everyday routine. Caregivers indicated the location type (e.g. home, kid’s
school, friends’ and relatives’ houses, work, grocery stores, etc.) and when
they typically spend time at the location (e.g. daytime, nighttime, weekdays,
weekends). They also answered multiple questions on each reported location to
assess the social climate. These location ratings are the key data used in
this current analysis to understand spatial variation in collective efficacy.
From the location assessments, we focus on three collective efficacy questions
posed as statements to which the caregiver rated agreement on a five-point
Likert scale: If someone was being threatened near [location], other people
around would come to their defense. You can trust people on the streets in the
area near [location]. There are usually people watching what’s happening in
the area near [location].111The first two statements had response options
”strongly agree,” ”agree,” ”neither agree nor disagree,” ”disagree,” and
”strongly disagree,” while the last statement had responses ”never,” ”almost
never,” ”sometimes,” ”fairly often” and ”very often.” Each response was
displayed in the order listed above next to the corresponding numeric value
1-5. In fitting the model we recoded the responses so that 1 corresponded to
the least affirmative value (never or strongly disagree) and 5 the most
affirmative value (very often and strongly agree) for all three components. We
refer to the responses to these questions and their corresponding latent
spatial processes, respectively, as ”defense,” ”trust,” and ”observation.”
While location ratings outside of the study area were provided, we restrict
our analysis to the 4526 locations within the I-270 belt loop. Out of the 1369
caregivers in the data set, the number of reports by a single caregiver ranged
from 1 to 27, with an average of 6.8 reports. With a few exceptions222Due to a
technical error in the administration of the survey, the question about trust
was omitted from the initial surveys resulting in fewer ratings for the trust
component as compared to defense and observation., all reports included
ratings on all three components of collective efficacy along with the time of
day the location is visited (daytime, nighttime, or both) and the day of week
(weekday, weekend, or both). Caregivers were asked to rate their home location
for both daytime and nighttime separately.

Figure 1: Survey responses for the defense (top row), trust (middle row), and
observation (bottom row) components plotted for residential locations (left
column) and nonresidential locations (right column). The gray lines show
census tract boundaries and the orange line is the 270 belt loop. Rated
locations have been jittered within census block group to preserve the
anonymity of study participants. Due to the jittering, some locations are
plotted outside the I-270 belt loop. The plots were created in R (R Core Team,
2022; Pebesma, 2018; Wickham et al., 2019).

Figure 1 shows the locations to which the ratings of defense, trust, and
observation correspond, plotted separately for residential and nonresidential
locations (as derived from the caregiver’s report on the location type). For
the defense responses, the residential ratings are consistent with
expectations if a spatially-dependent process underlies the ordinal ratings.
Spatial clustering of higher defense ratings in residential areas are apparent
in the northwest quadrant and a pocket of residences relative to the lower
ratings of defense elsewhere in the city. Nonresidential locations also appear
to exhibit the same general spatial patterning, yet not as well visually
pronounced. While it is not always appropriate to draw conclusions about the
strength of spatial dependence in a spatially-structured process underlying
spatially-resolved oridinal data (see Sections 3.1 and 3.3 for a discussion of
this issue), the visual differences are indicative of the process underlying
the ratings.

A similar pattern is evident in Figure 1 for the trust and observation
components – the spatial patterning of ratings for residential and
nonresidential locations are distinct. For trust (middle row), both
residential and nonresidential ratings exhibit clear spatial patterning, but
the pattern is not identical across both location types. In contrast, the
observation ratings (bottom row) visually show less of a distinctive spatial
pattern. It is difficult from visual inspection of these maps alone to draw
conclusions about the underlying spatial process that give rise to the
observed ratings across residential and nonresidential locations. Yet, for all
three components, it appears that allowing the mean process to vary as a
parametric function of observed, spatially-referenced covariates would not
capture the differences in the spatial patterning of observed ratings across
residential and nonresidential locations.

###  2.2 Land-use data

In order to account for the effects of land use on differences in the spatial
dependence structure of the components of collective efficacy across the study
region when predicting collective efficacy components at unrated locations, we
require the land-use designation of all prediction locations, which is not
available in the AHDC data. In particular, caregivers’ reported location types
are straightforward to classify as residential and nonresidential, however, we
need an objective and systematic method to categorize all other locations in
the study area not reported on. To this end, we utilized the May 2014 parcel
data for the study area, downloaded from the Franklin County Auditor’s
publicly accessible FTP site (, 2021). Each parcel has a designated tax code
which we classified into one of three categories: residential, nonresidential,
and other. We introduce the “other” category to distinguish nonresidential
locations with minimal social presence, such as warehouses, quarries, parking
lots and industrial centers, from locations where a social presence is
expected. Figure 2 shows the partition of the study area into residential,
nonresidential, and other categories as informed by the taxcodes of the parcel
data. Most locations in the AHDC data lie within the boundaries of a parcel
and received the same classification as the parcel within which they lie. The
coordinates of other AHDC-reported locations fall on the road network and do
not lie within the boundaries of a parcel. For these locations we assigned the
land-use category of the nearest parcel. After initial assignments were made,
we performed a second sweep of the data to verify that the reported caregiver
location type matched the assigned parcel land-use type. Incongruencies
existed between reported location types and assigned land-use types (e.g.
location reported by caregiver as “neighborhood,” but assigned to parcel with
tax designation “warehouse” and land-use category “other”). For such
locations, we reassigned the land-use category if any of the next 5 nearest
parcels had a tax code and land-use assignment that was coherent with the
reported location type. We deemed locations without a nearby congruent parcel
as a reporting error and did not include them in the analysis. Table 1 gives
the total number of ratings and locations used in our analysis broken down by
collective efficacy component and land-use category.

Table 1: Breakdown of the number of locations, $m$, and ratings, $n$, by collective efficacy component. The subscripts on $m$ and $n$ further breakdown location and ratings counts by land-use type, with 1 for residential, 2 for nonresidential and 3 for “other”. | Defense | Trust | Observation
---|---|---|---
$m$ | 3867 | 3040 | 3826
$m_{1}$ | 1770 | 1555 | 1777
$m_{2}$ | 2045 | 1447 | 1997
$m_{3}$ | 52 | 38 | 52
$n$ | 7580 | 5865 | 7548
$n_{1}$ | 2906 | 2653 | 2917
$n_{2}$ | 4616 | 3170 | 4573
$n_{3}$ | 58 | 42 | 58
Figure 2: Partition of the study region into land-use categories. Each cell
color is determined by the tax code of the nearest parcel

###  2.3 Preliminary exploratory modeling

To better understand the difference in spatial variation in collective
efficacy components by land-use categories, we fit separate univariate
Bayesian spatial ordinal regression models to each component of collective
efficacy by land-use category (excluding the “other” category which has only
approximately 50 ratings for each component). That is, we fit six models
corresponding to the three components of collective efficacy (defense, trust,
and observation) by the two land-use categories (residential and
nonresidential, as determined by the parcel data). Section 3 provides details
on the specification of spatial ordinal regression models but here we focus on
inferences on key parameters describing the nature of the spatial dependence
as a justification of our proposed model. Because Bayesian spatial ordinal
regression models are partially identifiable (discussed in Section 3.3), we
clarify that inferences on the spatial parameters cannot be interpreted as
characterizing the true latent continuously-indexed spatial process; rather,
these inferences serve as evidence of distinct spatial patterns by land-use
category. With this caveat in mind, let $\boldsymbol{\theta}=(\phi,\tau^{2})$
be the vector of spatial dependence parameters characterizing the exponential
covariance function of the latent spatial process, where $\phi$ represents the
spatial range and $\tau^{2}$ is the parameter that governs the proportion of
variance due to spatial variation, as opposed to independent error (in the
Bayesian ordinal regression model the variance of the independent error is
fixed at one). For all models, we adjust for the time-of-day (daytime,
nighttime, or both) and the day-of-week (weekday only, weekend only, or mixed)
variables in the mean function as indicated by the caregiver in the survey
response. We obtained draws from the posterior distribution using the MCMC
scheme described in Section 3.5.1.

Table 2: Estimates of the marginal posterior mean of the covariance parameters from Bayesian spatial ordinal models for the defense, trust, and observation components fitted separately to the residential and nonresidential ratings. | Defense |  | Trust |  | Observation
---|---|---|---|---|---
| $\phi$ | $\tau^{2}$ |  | $\phi$ | $\tau^{2}$ |  | $\phi$ | $\tau^{2}$
Residential | 334.29 | 2.08 |  | 443.68 | 3.03 |  | 2156.17 | 2.06
Nonresidential | 53.28 | 0.15 |  | 42.63 | 0.48 |  | 77.02 | 0.16

Table 2 shows the estimated mean of the marginal posterior distribution of the
spatial range parameter, $\phi$, and spatial proportion-of-variance parameter,
$\tau^{2}$, for all six fitted preliminary models. The parameter estimates for
the models fitted to the residential data differ from the estimates from the
nonresidential data across all three collective efficacy components. A higher
value of $\phi$ indicates a shorter range of spatial dependence and a high
value of $\tau^{2}$ indicates a larger proportion of the variance is
attributed to the spatial dependence rather than independent random error. We
plot the estimated posterior mean correlation function with 95% pointwise
credible intervals in Figure 3. For each pair of posterior samples of $\phi$
and $\tau^{2}$, we calculated the correlation function for a grid of distances
from 0 to 5 miles and plot the pointwise mean and lower/upper bounds of a 95%
credible interval. For the defense component, the estimated correlation
function for the residential data is characterized by shorter term spatial
dependence but a greater proportion of variation attributable to spatial
dependence. In contrast, the correlation function for the nonresidential data
of the defense component exhibits longer spatial dependence and a much larger
proportion of independent error. A similar pattern characterizes the
difference between the correlation functions for the trust and observation
components. These differences in the estimated correlation functions between
the two types of locations suggest that a model which accounts for the
difference in spatial dependence between locations of different land-use types
will better capture the overall spatial dependence structure of the data.

Figure 3: Estimated posterior mean correlation functions from Bayesian spatial
ordinal models for residential and nonresidential locations of the defense,
trust, and observation component, plotted respectively from left to right.
Dashed lines indicate lowers and upper bounds of 95% posterior point-wise
credible intervals.

The spatial dependence structure of the model has implications for the main
goal of our analysis: to predict levels of each component of collective
efficacy at unobserved locations. A spatial generalized linear mixed model
with a single latent stationary spatial process would most likely lead to a
spatial dependence parameter estimate somewhere between the distinct land-use
specific spatial dependence parameters shown in Table 2. We show this is
indeed the case in Section 4 for all three components of collective efficacy.
By expanding the dimension of the response to include a latent spatial process
for each land-use type, we can better separate the distinct spatial dependence
structures of each land-use process. This will allow us to have a better
understanding of the spatial process of each component of collective efficacy
at a finer resolution within census tract or block group. Furthermore, we can
smooth over land-use boundaries through the cross correlations between the
latent processes corresponding to land-use types.

##  3 Land-use filtering model

In this section, we introduce a novel land-use filtering model, a Bayesian
spatial generalized linear mixed model that incorporates dimension expansion
and multivariate spatial modeling strategies. We begin with a motivating
schematic to demonstrate how land-use filtering captures nonstationary
behavior in a latent univariate spatial process. Throughout this section we
omit notation references to the three different components of collective
efficacy since we will fit the same model to each component separately (see
Section 6 for more details on the rationale for fitting the components
separately).

###  3.1 Motivation

When the study area can be partitioned into different land-use categories (or
some other spatial partition), we can expand the dimension of the latent
spatial process in a spatial generalized linear mixed model to allow distinct
spatial dependence structures for each land-use category. To motivate the
land-use filtering model developed formally in the next section, consider the
following generative model illustrated in Figure 4. We generate from a mean-
zero latent multivariate spatial process with a known cross-covariance
function on a fine grid of locations distributed across a two-dimensional
region (a unit square in Figure 4). The assumed marginal spatial correlation
functions for this illustrative example are shown on the left panel. We then
discard the components of each sample that do not correspond to the land-use
category of a particular location, producing the top three plots in the center
column of Figure 4. The bottom plot shows the filtered process created by
combining the land-use specific component processes at each location. The
resulting image is fairly locally smooth across land-use boundaries due to the
cross-component dependence of the latent process, but the spatial dependence
structure is not identical across space. Lastly, we assign an ordinal rating
to a set of random locations within the study area based on the value of the
latent process at the locations. That is, these ratings are based on the non-
discarded component of the simulated random vector associated with that
location (top three plots on the right column). The simulated ratings (bottom
right) arising from the filtered latent process (bottom center) are analogous
to the AHDC collective efficacy ratings.

Figure 4: A demonstration of a generative model for ordinal data under land-
use filtering. See Section 3.1 for a discussion.

Note that in Figure 4 when the ratings across all three land-use categories
are plotted together (bottom right column), it is difficult to visualize
different spatial dependence structures which are (more) distinguishable when
the latent processes are plotted together. By expanding the univariate
response to allow for distinct spatial correlation functions for each land-use
category and dependence across land-use categories, we can improve model fit
for the AHDC collective efficacy ratings.

###  3.2 Land-use filtering ordinal regression

In specifying our land-use filtering model, we consider only a single
component of collective efficacy. Each of the AHDC caregivers provides their
rating on a five-point ordinal scale, where a rating of $K\equiv 5$ is the
“best” and 1 is the “worst.” To facilitate specification of our model, we
introduce notation for the ratings, with the primary index indicating the
location of the rating. We let $Y_{ij}$ denote the $j$th rating of the $i$th
location, where $\mathbf{s}_{i}$ is a unique location in the study area,
$\mathcal{S}$, for $i=1,\dots,m$ and $j=1,\dots,n_{i}$. The constant $m$ is
the number of unique locations, and $n_{i}$ is the total number of caregivers
who rated the $i$th location. We denote the total number of ratings as
$n=\sum_{i=1}^{m}n_{i}$, and the collection of all ratings as
$\mathcal{Y}=\\{Y_{ij}:i=1,\dots,m;j=1,\dots n_{i}\\}$.

We model the elements of $\mathcal{Y}$ using a Bayesian ordinal probit
regression model, specified using the well known data augmentation scheme of
Albert and Chib (1993, 1997), which has been extended to the spatial setting
by De Oliveira (1997, 2000); Higgs and Hoeting (2010); Schliep and Hoeting
(2015); Berrett and Calder (2012, 2016). Under this scheme, we introduce
continuous latent variables $\mathrm{Z}_{ij}^{*}$ defined at each of the $m$
locations and for each $n_{i}$ ratings and unknown cut points,
$\gamma_{1},\dots,\gamma_{K-1}$ such that

$Y_{ij}=\begin{cases}K&\mbox{if }\gamma_{K-1}<\mathrm{Z}^{*}_{ij}\\\
k&\mbox{if }\gamma_{k-1}<\mathrm{Z}^{*}_{ij}<\gamma_{k},\mbox{ for
}k=2,\dots,K-1\\\ 1&\mbox{if }\mathrm{Z}^{*}_{ij}<\gamma_{1}\end{cases}.$

Unlike the setting discussed in Higgs and Hoeting (2010), we have multiple
latent random variables defined at each location due to the fact that some
locations are rated by multiple individuals. Each of these random variables is
related to a single latent spatial process $\tilde{Z}(\mathbf{s})$ defined for
all $\mathbf{s}\in\mathcal{S}.$ We assume that the $\mathrm{Z}_{ij}^{*}$s are
conditionally independent given the latent spatial process
$\tilde{Z}(\mathbf{s}_{i})$, and parameters $\boldsymbol{\beta}$,
$\boldsymbol{\theta}$, and $\boldsymbol{\sigma}^{2}$:

$\mathrm{Z}^{*}_{ij}|\tilde{Z}(\mathbf{s}_{i}),\boldsymbol{\beta},\boldsymbol{\theta},\tilde{\boldsymbol{\sigma}}^{2}\sim\mathrm{N}(\mathbf{x}_{ij}^{\prime}\boldsymbol{\beta}+\tilde{Z}(\mathbf{s}_{i}),\;\tilde{\sigma}^{2}_{g(\mathbf{s}_{i})}).\\\
$

Here, $\mathbf{x}_{ij}$ is a $p\times 1$ vector of covariates associated with
the $j$th rating of location $i$ and $\boldsymbol{\beta}$ is the corresponding
$p\times 1$ vector of regression coefficients. The vector
$\tilde{\boldsymbol{\sigma}}^{2}$ contains the variance parameters for the
independent error of each land-use process and
$g(\mathbf{s}):\mathcal{S}\rightarrow\\{1,\dots,Q\\}$ is a function that
returns the index of the land-use category of location $\mathbf{s}$. We use
the tilde notation to highlight the fact that this latent process is
constrained so that the total error variance of the $\mathrm{Z}_{ij}^{*}$s is
equal to one. Constraining the latent spatial process in ordinal regression
model can aid in the identifiability of the covariance parameters
$\boldsymbol{\theta}$ (Schliep and Hoeting, 2015; Yan et al., 2007). We
discuss additional identification considerations for spatial ordinal
regression models in Section 3.3. The vector $\boldsymbol{\theta}$ denotes the
covariance parameters associated with the spatial process
$\tilde{Z}(\mathbf{s})$, on which $\tilde{\boldsymbol{\sigma}}^{2}$ also
depends, due to the constraint on the total error variance of
$\mathrm{Z}_{ij}^{*}$.

Before defining the constrained $\tilde{Z}(\mathbf{s})$ and a model for its
unconstrained analogue $Z(\mathbf{s})$, we introduce a $Q$-dimensional spatial
process $\boldsymbol{\eta}(\mathbf{s})$, where $Q$ is the number of land-use
categories in the spatial partition. We define $\boldsymbol{\eta}(\mathbf{s})$
to be a multivariate Gaussian process with mean zero and parametric covariance
function
$\boldsymbol{\Sigma}^{(\boldsymbol{\eta})}_{\mathbf{s},\mathbf{s}^{\prime}}(\boldsymbol{\theta})$
defined by the linear model of coregionalization (LMC). The LMC remains the
standard approach to create covariance functions in multivariate spatial
modeling (Goulard and Voltz, 1992; Grzebyk and Wackernagel, 1994; Schmidt and
Gelfand, 2003; Wackernagel, 2003; Gelfand et al., 2004). While alternatives to
the LMC exist (e.g., Gneiting, Kleiber and Schlather (2010) and Apanasovich,
Genton and Sun (2012) propose Matérn cross-covariance functions for
multivariate random fields), De Oliveira (1997, 2000) has shown that
smoothness parameters (such as the smoothness parameter in the Matérn
covariance function) are near nonidentifiable in spatial ordinal regression
models. Since our $\boldsymbol{\eta}$ will be a component of a spatial ordinal
regression model, we prefer the LMC specification.

The LMC is a constructive approach for modeling a multivariate spatial process
as a linear combination of independent spatial processes. Let
$\mathbf{w}(\mathbf{s})$ be a $Q$-dimensional spatial process comprised of
independent component processes
$\mathbf{w}_{1},\dots,\mathbf{w}_{Q}(\mathbf{s})$. Each
$\mathbf{w}_{q}(\mathbf{s})$ is a Gaussian process with mean 0, variance 1,
and correlation function $\rho_{q}(\mathbf{s},\mathbf{s}^{\prime};\phi_{q})$,
where $\phi_{q}$ are unknown parameters. Again, since smoothness parameters
are near nonidentifiable, we specify
$\rho_{q}(\mathbf{s},\mathbf{s}^{\prime};\phi_{q})$ as an exponential
correlation function,

$\rho_{q}(\mathbf{s},\mathbf{s}^{\prime},\phi_{q})=\exp(-\phi_{q}||\mathbf{s}-\mathbf{s}^{\prime}||),$

but note that alternate correlation functions could be used instead. In the
LMC specification, we assume that
$\boldsymbol{\eta}(\mathbf{s})=\mathbf{A}\mathbf{w}(\mathbf{s})$, where
$\mathbf{A}$ is a $Q\times Q$ lower triangular matrix with $Q\times(Q-1)/2$
unknown parameters. It follows directly that the cross covariance for
$\boldsymbol{\eta}(\mathbf{s})$ and $\boldsymbol{\eta}(\mathbf{s}^{\prime})$
is

$\boldsymbol{\Sigma}^{(\boldsymbol{\eta})}_{\mathbf{s},\mathbf{s}^{\prime}}(\boldsymbol{\theta})=\mathbf{A}\boldsymbol{\Gamma}(\mathbf{s},\mathbf{s}^{\prime})\mathbf{A}^{\prime},$

where $\boldsymbol{\theta}$ consists of the elements of $\boldsymbol{\phi}$
and unknown elements of $\mathbf{A}$ and
$\boldsymbol{\Gamma}(\mathbf{s},\mathbf{s}^{\prime})$ is a matrix valued
function with $\rho_{q}(\mathbf{s},\mathbf{s}^{\prime},\phi_{q})$ on the $q$th
diagonal.

We now relate the Q-dimensional process $\boldsymbol{\eta}(\mathbf{s})$ to the
univariate process $Z(\mathbf{s})$ using land-use filtering (i.e., extracting
the component of $\boldsymbol{\eta}(\mathbf{s})$ corresponding to the land-use
category of $\mathbf{s}$). Formally, we define

$Z(\mathbf{s})=f(\boldsymbol{\eta}(\mathbf{s}),g(\mathbf{s}))=\eta_{g(\mathbf{s})}(\mathbf{s})+\delta_{g(\mathbf{s})},$

where $\boldsymbol{\delta}=(\delta_{1},\dots,\delta_{Q})$ is a vector of
unknown mean shift parameters defining the mean of the $Z(\mathbf{s})$
process. Writing
$Z(\mathbf{s})=\mathbf{a}_{g(\mathbf{s})}^{\prime}\mathbf{w}(\mathbf{s})+\delta_{g(\mathbf{s})}$,
where $\mathbf{a}_{q}$ is the $q$th row of $\mathbf{A}$, it is clear that
$Z(\mathbf{s})$ is a Gaussian process with mean function $\delta(\mathbf{s})$,
where

$\delta(\mathbf{s})=\delta_{g(\mathbf{s})},$

and cross-covariance function

$\boldsymbol{\Sigma}^{(Z)}_{\mathbf{s},\mathbf{s}^{\prime}}(\boldsymbol{\theta})=\mathbf{a}_{g(\mathbf{s})}^{\prime}\boldsymbol{\Gamma}(\mathbf{s},\mathbf{s}^{\prime})\mathbf{a}_{g(\mathbf{s}^{\prime})}.$

between any two locations $\mathbf{s}$ and $\mathbf{s}^{\prime}$ with land-use
types $g(\mathbf{s})$ and $g(\mathbf{s}^{\prime})$. Note that $Z(\mathbf{s})$
is nonstationary with both a mean function $\boldsymbol{\delta}(\mathbf{s})$
and covariance function
$\boldsymbol{\Sigma}^{(Z)}_{\mathbf{s},\mathbf{s}^{\prime}}(\boldsymbol{\theta})$
that vary across land-use categories.

To complete the specification of the model, we need to define the covariance
function for the constrained $\tilde{Z}(\mathbf{s})$, which allows the
$Z^{*}_{ij}$s to have total variance of one. First we note that the variance
of the unconstrained, land-use-category-specific process
$\eta_{q}(\mathbf{s})+\delta_{q}$ is $||\mathbf{a}_{q}||^{2}$. It follows then
that the constrained covariance is

$\tilde{\boldsymbol{\Sigma}}^{(\tilde{Z})}_{\mathbf{s},\mathbf{s}^{\prime}}(\boldsymbol{\theta}^{*})=\frac{\mathbf{a}_{g(\mathbf{s})}^{\prime}\boldsymbol{\Gamma}(\mathbf{s},\mathbf{s}^{\prime})\mathbf{a}_{g(\mathbf{s}^{\prime})}}{\sqrt{||\mathbf{a}_{g(\mathbf{s})}||^{2}+\sigma^{2}_{g(\mathbf{s})}}\sqrt{||\mathbf{a}_{g(\mathbf{s}^{\prime})}||^{2}+\sigma^{2}_{g(\mathbf{s}^{\prime})}}},$
(1)

where
$\boldsymbol{\theta}^{*}=\\{\mathrm{vec}(\mathbf{A}),\boldsymbol{\phi},\boldsymbol{\sigma}^{2}\\}$.
Additionally, the independent errors are also constrained so that

$\tilde{\sigma}^{2}_{g(\mathbf{s})}=\frac{\sigma^{2}_{g(\mathbf{s})}}{||\mathbf{a}_{g(\mathbf{s})}||^{2}+\sigma^{2}_{g(\mathbf{s})}}.$

The resulting constrained process, $\tilde{Z}(\mathbf{s})$, is a nonstationary
Gaussian process

$\tilde{Z}(\mathbf{s})|\delta(\mathbf{s}),\boldsymbol{\theta}^{*}\sim\mathrm{GP}(\delta(\mathbf{s}),\tilde{\boldsymbol{\Sigma}}^{(\tilde{Z})}_{\mathbf{s},\mathbf{s}^{\prime}}(\boldsymbol{\theta}^{*})),$

with mean function, $\boldsymbol{\delta}(\mathbf{s})$, and covariance
function,
$\tilde{\boldsymbol{\Sigma}}^{(\tilde{Z})}_{\mathbf{s},\mathbf{s}^{\prime}}(\boldsymbol{\theta}^{*})$.

###  3.3 Identifiability of model parameters

It is well known that Bayesian probit regression models for ordinal and binary
data (and subsequently the spatial extensions of these models) contain
nonidentifiable parameters. That is, there is not a one-to-one mapping between
the parameters and the value of the likelihood function. We describe below the
necessary constraints to obtain likelihood identifiability in the land-use
filtering model.

First, a global intercept (or cell means intercepts) and the cut points
$\boldsymbol{\gamma}$ cannot be jointly identified. To address this
identification issue we fix $\delta_{1}=0$, the mean of the first land-use
process, and estimate the remaining $\delta_{q}$ for $q>1$. Additionally, we
exclude a global intercept from $\boldsymbol{\beta}$, the vector of regression
coefficients. Alternatively, one could fix one of the cut points, for example
$\gamma_{1}=0$, which would allow the estimation of a global intercept in
$\boldsymbol{\beta}$ or the estimation of $\delta_{1}$, but not both. As a
last consideration on the cut points, in a standard multivariate probit model,
one could use a different set of cut points for each multivariate process,
although this would preclude the use of land-use specific intercepts. We opt
for a single set of cut points across the land-use processes so that
interpretation of the latent variable $\tilde{Z}(\mathbf{s})$ is consistent
across land-use categories.

In Bayesian spatial probit regression models with an independent error term,
the variance of the independent error and the regression coefficients
$\boldsymbol{\beta}$ are identified up to a multiplicative constant. In the
hierarchical model specified above, one could drop the constraint on the
spatial random effects $\mathbf{Z}(\mathbf{s})$ and fix each $\sigma^{2}_{q}$
to obtain likelihood identifiability. We instead choose to constrain the total
error variance of $\mathrm{Z}^{*}$ to be one, as this can aid the
identifiability of the covariance parameters (Schliep and Hoeting, 2015) and
facilitates a block update of the covariance parameters in our MCMC algorithm
as described in Section 3.5.1. Under this scaling regime of equation 1, we are
limited to estimating the proportion of the variance due to spatial
dependence, and again the $\sigma_{q}^{2}$’s must all be fixed. By fixing
$\sigma^{2}_{q}=1$ for all $q$, we can express the proportion of variance due
to spatial correlation for each land-use category $q$ as
$\frac{||\mathbf{a}_{q}||^{2}}{||\mathbf{a}_{q}||^{2}+1}$.

As the caregiver dataset has replication at the location level, we are
required to include independent error for the latent process of the ratings.
For a setting without replication, one could drop the independent error
assumption and let

$\mathrm{Z}^{*}_{i}=\mathbf{x}_{i}^{\prime}\boldsymbol{\beta}+\tilde{Z}(\mathbf{s}_{i}).$

Now we constrain the variance of $\tilde{Z}(\mathbf{s})$ to be one, as again
an unconstrained covariance and $\boldsymbol{\beta}$ are nonidentified. The
constrained covariance for the model without independent error results in the
correlation function

$\boldsymbol{R}^{(\tilde{Z})}_{\mathbf{s},\mathbf{s}^{\prime}}(\boldsymbol{\theta})=\frac{\mathbf{a}_{g(\mathbf{s})}^{\prime}\boldsymbol{\Gamma}(\mathbf{s},\mathbf{s}^{\prime})\mathbf{a}_{g(\mathbf{s}^{\prime})}}{||\mathbf{a}_{g(\mathbf{s})}||\,||\mathbf{a}_{g(\mathbf{s}^{\prime})}||}.$

Lastly, we discuss another identifiability concern that is uniquely Bayesian,
that of partial identifiability. Partially identifiable models are
characterized by posterior distributions that do not converge to a point mass
as the sample size increases to infinity, yet are still informed by the data
and differ from the prior distribution (Gustafson, 2015). Partially identified
models contain partially or weakly identified parameters that are typified by
substantive flat regions in the posterior or by a posterior that critically
depends on the prior (Li, Ding and Mealli, 2022). Identifiability, in this
context, is on a continuum reflecting the strength of learning from the data.
Our model is partially identifiable in that our observed data is a corrupted
version of the collective efficacy process that we are trying to predict. In
particular, the covariance parameters are weakly identifiable. As a
consequence, we cannot expect to precisely learn about these parameters even
when the data are simulated from the model and the sample size is large.
Section 5 investigates the implication of the partial identifiability of our
model in terms of out-of-sample predictive performance.

###  3.4 Implications of the LMC specification

We now explore the implications of the LMC specification on the cross-
covariance function of the multivariate latent process. In the continuous
response setting, Gneiting, Kleiber and Schlather (2010) note that when all
elements of the matrix $\mathbf{A}$ are treated as unknown parameters (i.e. no
structural zeros), the smoothness of each latent process is determined by the
roughest latent component and that by fixing the upper right elements of the
matrix $\mathbf{A}$ to be zero, distinct smoothness properties can be
estimated. In a similar vein, we find that the structural zeros also help in
identifying distinct spatial range parameters for each component process.

Imposing this structure of $\mathbf{A}$, however, does indirectly place
restrictions on the cross-covariance function. Consider the setting where the
dimension of the latent process, $Q$, is three and $\mathbf{A}$ is lower
triangular such that $a_{qq^{\prime}}=0$ for all $q^{\prime}>q$, where
$a_{qq^{\prime}}$ denotes the $(q,q^{\prime})$ element of $\mathbf{A}$. Let
$\tilde{\boldsymbol{\eta}}(\mathbf{s})$ be the constrained multivariate
process defined by the LMC with the constraint given in Equation 1. The cross-
covariance matrix-valued function for $\tilde{\boldsymbol{\eta}}(\mathbf{s})$
and $\tilde{\boldsymbol{\eta}}(\mathbf{s}^{\prime})$ can be written as

$\begin{bmatrix}\frac{a_{11}^{2}\rho_{1}}{||\mathbf{a}_{1}||^{2}+\sigma^{2}_{1}}&\frac{a_{11}a_{21}\rho_{1}}{\sqrt{||\mathbf{a}_{1}||^{2}+\sigma^{2}_{1}}\sqrt{||\mathbf{a}_{2}||^{2}+\sigma^{2}_{2}}}&\frac{a_{11}a_{31}\rho_{1}}{\sqrt{||\mathbf{a}_{1}||^{2}+\sigma^{2}_{1}}\sqrt{||\mathbf{a}_{3}||^{2}+\sigma^{2}_{3}}}\\\
\frac{a_{11}a_{21}\rho_{1}}{\sqrt{||\mathbf{a}_{1}||^{2}+\sigma^{2}_{1}}\sqrt{||\mathbf{a}_{2}||^{2}+\sigma^{2}_{2}}}&\frac{a_{21}^{2}\rho_{1}+a_{22}^{2}\rho_{2}}{||\mathbf{a}_{2}||^{2}+\sigma^{2}_{2}}&\frac{a_{21}a_{31}\rho_{1}+a_{22}a_{32}\rho_{2}}{\sqrt{||\mathbf{a}_{1}||^{2}+\sigma^{2}_{1}}\sqrt{||\mathbf{a}_{3}||^{2}+\sigma^{2}_{3}}}\\\
\frac{a_{11}a_{31}\rho_{1}}{\sqrt{||\mathbf{a}_{1}||^{2}+\sigma^{2}_{1}}\sqrt{||\mathbf{a}_{3}||^{2}+\sigma^{2}_{3}}}&\frac{a_{21}a_{31}\rho_{1}+a_{22}a_{32}\rho_{2}}{\sqrt{||\mathbf{a}_{1}||^{2}+\sigma^{2}_{1}}\sqrt{||\mathbf{a}_{3}||^{2}+\sigma^{2}_{3}}}&\frac{a_{31}^{2}\rho_{1}+a_{32}^{2}\rho_{2}+a_{33}^{2}\rho_{3}}{||\mathbf{a}_{3}||^{2}+\sigma^{2}_{3}}\end{bmatrix},$
(2)

where $\rho_{q}=\rho(\mathbf{s},\mathbf{s}^{\prime},\phi_{q})$. From this
expression for the cross-correlation function in Equation 2, it is clear that
the cross-covariance between the first latent process and the other two only
depends on $\rho_{1}$. That is, the structural zeros in $\mathbf{A}$ restrict
the cross-covariance function between the first and subsequent latent process
to only depend on the spatial dependence structure of the first spatial
process. It immediately follows that the lower-triangular form of $\mathbf{A}$
implies that the ordering of the land-use categories will affect the estimated
cross-covariance function. Additionally, we note that ordering of the
components can affect the (partial) identifiability of distinct spatial range
parameters of the spatial processes. We discuss our strategy for ordering
land-use categories in Section 4.

Additionally, the implied cross-covariance function from the LMC model also
has implications for the trade off between independent error and spatial
variance when all $\sigma_{q}^{2}$s are set to one. For each land-use process
the proportion of independent error is $\frac{1}{||\mathbf{a}_{q}||^{2}+1}$. A
land-use process with greater independent error will have proportionally less
spatial error. Because the total error variance is constrained to one, the
land-use processes $\eta_{q}(\mathbf{s})+\delta_{q}$ will have smaller spatial
variance for land-use types with a greater proportion of independent error.

###  3.5 Model fitting

The full Bayesian model is specified by placing priors on the covariance
parameters $\boldsymbol{\theta}^{*}$, and mean function parameters
$\boldsymbol{\delta}$ and $\boldsymbol{\beta}$. A normal prior for
$\boldsymbol{\delta}$ and $\boldsymbol{\beta}$ is conditionally conjugate in
the data augmented probit model. The diagonal elements of $\mathbf{A}$,
$\boldsymbol{\phi}$, and $\boldsymbol{\sigma}^{2}$ have support on the
positive real line, while the off diagonal elements of $\mathbf{A}$ have
support on the real line.

We wrote a custom MCMC algorthim in Julia to fit the model (Bezanson et al.,
2017; Rackauckas and Nie, 2017; Lin et al., 2019; Besançon et al., 2021; Frigo
and Johnson, 2005; Mogensen and Riseth, 2018). As MCMC is computationally time
consuming and burdensome we also propose an approximation method that does not
utilize the data augmentation framework, but rather treats the ordinal ratings
as a Gaussian outcome. We fit the approximate model by direct maximization of
the approximate posterior. We say “approximate” as it is not appropriate to
assume ordinal data follow a normal distribution, but for the AHDC collective
efficacy data spatial predictions from this approximate model are nearly
identical to those produced by the land-use filtering model and can be
obtained much faster (computation time is hours instead of days).

To facilitate the discussion of each method used to fit the model we introduce
vector notation. Let $\mathbf{Y}$ be a $n\times 1$ vector of ratings and
$\mathbf{Z}^{*}$ be the $n\times 1$ latent response vector with $Y_{r}$ and
$\mathrm{Z}^{*}_{r}$ denoting the $r$th rating and latent variable. The vector
$\tilde{\mathbf{Z}}=(\tilde{\mathrm{Z}}_{1},\dots,\tilde{\mathrm{Z}}_{m})^{\prime}$
contains the spatial random effects for each location, where
$\tilde{\mathrm{Z}}_{i}=\tilde{Z}(\mathbf{s}_{i})$. The $r$th row of the
design matrix, $\mathbf{X}$, contains the covariates $\mathbf{x}_{r}$ of the
$r$th rating. $\mathbf{H}$ is a $n\times m$ matrix of ones and zeros that
associates the spatial random effects to the corresponding locations of each
rating. The $n\times 1$ vector $\boldsymbol{\epsilon}$ of independent errors
have a unique constrained variance, $\tilde{\sigma}_{q}^{2}$, for each land-
use type $q$. We then write the model as

$\mathbf{Z}^{*}=\mathbf{X}\boldsymbol{\beta}+\mathbf{H}\tilde{\mathbf{Z}}+\tilde{\boldsymbol{\epsilon}},$

where vector of spatial random effects, $\tilde{\mathbf{Z}}$, has a
multivariate normal distribution,

$\tilde{\mathbf{Z}}|\boldsymbol{\delta},\boldsymbol{\theta}^{*}\sim\mathrm{MVN}(\mathbf{M}\boldsymbol{\delta},\tilde{\boldsymbol{\Sigma}}(\boldsymbol{\theta}^{*})),$

and $\mathbf{M}$ is a $m\times Q$ matrix of zeros and ones that associates the
$i$th location with the mean shift $\delta_{q}$ of the corresponding land-use
type. Lastly, the $\tilde{\boldsymbol{\epsilon}}$ are independent,

$\tilde{\boldsymbol{\epsilon}}|\boldsymbol{\theta}^{*}\sim\mathrm{N}(\mathbf{0},\tilde{\mathbf{D}}),$

with $\tilde{\sigma}_{q}^{2}$ on the $r$th diagonal of $\tilde{\mathbf{D}}$
for rating $r$ with land-use type $q$.

####  3.5.1 MCMC

To facilitate mixing of the MCMC algorithm, we integrate out the spatial
random effects. After doing so, we write the model as

$\mathbf{Z}^{*}=\mathbf{X}\boldsymbol{\beta}+\mathbf{M}^{*}\boldsymbol{\delta}+\boldsymbol{\nu}\\\
,$

with the error term $\boldsymbol{\nu}$ having a multivariate normal
distribution,

$\boldsymbol{\nu}|\boldsymbol{\theta}^{*}\sim\mathrm{N}(0,\boldsymbol{R}(\boldsymbol{\theta}^{*})),$

where $\mathbf{M}^{*}=\mathbf{H}\mathbf{M}$ and

$\boldsymbol{R}(\boldsymbol{\theta}^{*})=\mathbf{H}\tilde{\boldsymbol{\Sigma}}(\boldsymbol{\theta}^{*})\mathbf{H}^{\prime}+\tilde{\mathbf{D}}$
(3)

is a correlation matrix. First, we update the cut points and latent vector
$\mathbf{Z}^{*}$. The cut points are drawn from a uniform distribution
$\gamma_{k}\sim\mathrm{U}(\max(\mathrm{Z}_{r}^{*}|Y_{r}=k),\min(\mathrm{Z}_{r}^{*}|Y_{r}>k))$.
While alternative schemes have been proposed for updating the cut points,
these strategies prove to be impractical computationally for our relatively
large data set in the spatial setting (Albert and Chib, 1997; Cowles, 1996;
Higgs and Hoeting, 2010). Each $\mathrm{Z}_{r}^{*}$ in the latent vector
$\mathbf{Z}^{*}$ is drawn from a normal distribution obtained by conditioning
on $\boldsymbol{\beta}$, $\boldsymbol{\delta}$ and all other observations
$\mathbf{Z}^{*}_{-r}$ and truncated to support
$[\gamma_{Y_{r-1}},\gamma_{Y_{r}}]$. As the precision matrix,
$\boldsymbol{\Lambda}=\boldsymbol{R}(\boldsymbol{\theta}^{*})^{-1}$, is used
to evaluate the likelihood of $\mathbf{Z}^{*}$ in the posterior for a
Metropolis-Hastings update of $\boldsymbol{\theta}$, we also write the
variance of $\mathrm{Z}_{r}^{*}|\mathbf{Z}_{-r}^{*}$ in terms of the precision
to avoid additional inverse calculations:
$\mathrm{var}(\mathrm{Z}_{r}^{*}|\mathbf{Z}^{*}_{-r})=1/\boldsymbol{\Lambda}_{rr}$.
After obtaining a draw of $\mathbf{Z}^{*}$, the regression coefficients
$\boldsymbol{\beta}$ and mean shift $\boldsymbol{\delta}$ can be drawn jointly
from their full conditional distributions with Gibbs updates. We update
$\boldsymbol{\theta}$ jointly with a Metropolis-Hastings step evaluating the
posterior conditional on the current draw of $\mathbf{Z}^{*}$. We use a
multivariate normal proposal distribution for $\boldsymbol{\theta}$ centered
around the current draw of $\boldsymbol{\theta}$. We tune the covariance
matrix for the proposal by running the algorithm to obtain at least ten
thousand draws with a diagonal covariance. From this initial run, we scale the
covariance of the draws of $\boldsymbol{\theta}$ by $0.25$ and use the
resulting matrix as the covariance of our proposal distribution for the final
run of the MCMC.

####  3.5.2 Approximate maximum a posteriori estimation

While running the full MCMC algorithm allows us to assess uncertainty in the
model parameters, the algorithm is very time consuming to run due to the
matrix inversions required to compute the likelihood of the latent vector
$\mathbf{Z}^{*}$. As an alternative, we consider a strategy to estimate the
model parameters that finds the maximum a posteriori estimates of parameters
from a model that approximates the full Bayesian ordinal specification. The
approximate model treats the outcomes as Gaussian, so that

$\displaystyle\mathbf{Y}^{*}$
$\displaystyle=\mathbf{X}\boldsymbol{\beta}+\mathbf{M}^{*}\boldsymbol{\delta}+\boldsymbol{\nu}$
$\displaystyle\boldsymbol{\nu}|\boldsymbol{\theta}^{*}$
$\displaystyle\sim\mathrm{N}(\mathbf{0},\boldsymbol{\Sigma}(\boldsymbol{\theta}^{*})),$

where
$\boldsymbol{\Sigma}(\boldsymbol{\theta}^{*})=\mathbf{H}\boldsymbol{\Sigma}(\boldsymbol{\theta})\mathbf{H}^{\prime}+\mathbf{D}$.
The main difference between the ordinal model and the approximate model is
that the covariance of $\boldsymbol{\nu}$ is unconstrained. To simplify the
calculations required to find the a posteriori estimates under this
approximate model, we integrate out $\boldsymbol{\beta}$ and
$\boldsymbol{\delta}$ and only estimate the covariance parameters
$\boldsymbol{\theta}^{*}$. The marginal distribution of $\mathbf{Y}$ is then
given by

$\mathbf{Y}|\boldsymbol{\theta}^{*}\sim\mathrm{N}(\mathbf{0},\mathbf{X}\boldsymbol{\Sigma}_{\boldsymbol{\beta}}\mathbf{X}^{\prime}+\mathbf{M}^{*}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}{\mathbf{M}^{*}}^{\prime}+\mathbf{H}\boldsymbol{\Sigma}(\boldsymbol{\theta})\mathbf{H}^{\prime}+\mathbf{D}),$

where $\boldsymbol{\Sigma}_{\boldsymbol{\delta}}$ and
$\boldsymbol{\Sigma}_{\boldsymbol{\beta}}$ are the covariance from the priors
placed on $\boldsymbol{\delta}$ and $\boldsymbol{\beta}$. We optimize for
$\boldsymbol{\theta}^{*}$ by placing priors on $\boldsymbol{\theta}^{*}$ and
then finding the maximum a posteriori estimate of the approximate marginal
posterior.

###  3.6 Prediction of the latent collective efficacy process

While the land-use filtering model is built on the assumption that we only
observe a single land-use process at each location, when making predictions we
can predict the latent level of our response for all land-use categories at
any location. These predictions may be useful in visualizing how the spatial
processes differ from one another, however, predictions of land-use processes
at locations that do not match the true land-use category may not be
meaningful depending on the application at hand. It would not be correct, for
our data example, to predict the change in the latent level of defense if a
location changed land-use categories from residential to nonresidential.

We now detail how to obtain samples from the posterior predictive distribution
for the spatial process $\tilde{Z}(\mathbf{s})$ at any location. Let
$\bar{\mathcal{S}}$ be the ordered set of new locations at which we want to
obtain predictions of the spatial process $\tilde{Z}(\mathbf{s})$ for
$\mathbf{s}\in\bar{\mathcal{S}}$ for a component of collective efficacy. Let
$\bar{\mathbf{M}}$ be the matrix which associates prediction locations to a
land-use category. Additionally, let draws from the posterior distribution be
indexed by $b$. To obtain a draw $\tilde{\mathbf{Z}}^{[b]}$ from the posterior
predictive distribution we need to evaluate the following covariance matrices
plugging in the posterior samples of model parameters (denoted by superscript
$[b])$:
$\tilde{\boldsymbol{\Sigma}}_{\bar{\mathcal{S}},\mathcal{S}}(\boldsymbol{\theta}^{*[b]})$
and
$\tilde{\boldsymbol{\Sigma}}_{\bar{\mathcal{S}},\bar{\mathcal{S}}}(\boldsymbol{\theta}^{*[b]})$
from Equation 1 and
$\boldsymbol{R}_{\mathcal{S},\mathcal{S}}(\boldsymbol{\theta}^{*[b]})$ from
Equation 3. Then using properties of the multivariate normal distribution, we
draw $\tilde{\mathbf{Z}}^{[b]}$ from $\tilde{\mathbf{Z}}|\mathbf{Z}^{*[b]}$
from $\mathrm{N}(\boldsymbol{\mu}^{[b]},\boldsymbol{\Sigma}^{[b]})$, where

$\displaystyle\boldsymbol{\mu}^{[b]}$
$\displaystyle=\bar{\mathbf{M}}\boldsymbol{\delta}^{[b]}+\tilde{\boldsymbol{\Sigma}}_{\bar{\mathcal{S}},\mathcal{S}}(\boldsymbol{\theta}^{*[b]})\boldsymbol{R}_{\mathcal{S},\mathcal{S}}^{-1}(\boldsymbol{\theta}^{*[b]})(\mathbf{Z}^{*[b]}-\mathbf{X}\boldsymbol{\beta}^{[b]}-\mathbf{M}^{*}\boldsymbol{\delta}^{[b]})$
(4) $\displaystyle\boldsymbol{\Sigma}^{[b]}$
$\displaystyle=\tilde{\boldsymbol{\Sigma}}_{\bar{\mathcal{S}},\bar{\mathcal{S}}}(\boldsymbol{\theta}^{*[b]})-\tilde{\boldsymbol{\Sigma}}_{\bar{\mathcal{S}},\mathcal{S}}(\boldsymbol{\theta}^{*[b]})\boldsymbol{R}_{\mathcal{S},\mathcal{S}}^{-1}(\boldsymbol{\theta}^{*[b]})\tilde{\boldsymbol{\Sigma}}_{\mathcal{S},\bar{\mathcal{S}}}(\boldsymbol{\theta}^{*[b]}).$

The $\tilde{\mathbf{Z}}^{[b]}$s are samples from the posterior predictive
distribution $\tilde{\mathbf{Z}}|\mathcal{Y}$.

##  4 Analysis of the AHDC caregiver data

We fit our proposed land-use filtering ordinal regression model to the
caregiver ratings, separately for each of the three components of collective
efficacy. As noted in Section 3.4, the ordering of the land-use processes in
the LMC specification is consequential for identifiability of distinct spatial
range parameters and has implications on the cross-covariance between land-use
types. From our preliminary analyses in Section 2.3, we found that the
residential processes had shorter spatial range than the nonresidential latent
processes. Therefore, we selected residential as the first component in the
LMC specification of the multivariate process $\boldsymbol{\eta}(\mathbf{s})$
and nonresidential as the second. We then set the “other” category third, as
the “other” category has very few observations ($\approx 50$) for all three
components of collective efficacy. Moreover, we set the third latent component
process in the LMC specification for the covariance function to be spatially
independent (i.e. setting $\phi_{3}$ to positive infinity) as we found that
the range parameter of the third latent process, $\phi_{3}$, is weakly
identifiable for each component and MCMC chains did not appear to converge
when we attempted to estimate it. As we expect the components of collective
efficacy to differ in terms of the nature of the spatial dependence structure
between nonresidential and “other” land-use areas and the “other” category
covers a nonnegligible portion of the city, we could not justify collapsing
the nonresidential and “other” categories into a single land-use category.

In the design matrix, $\mathbf{X}$, we include time-of-day indicators
(daytime, nighttime, or mixed) and the day-of-week indicators (weekday only,
weekend only, or mixed) as reported by the survey participants. We allow the
fixed effects corresponding to time-of-day and day-of-the-week to differ by
land-use category. The baseline category is residential locations rated during
the daytime. We do not include a day-of-week effect for residential locations
as nearly all reports were day-of-week mixed. For the “other” category, we
collapse day-of–week categories weekend-only and mixed and time-of-day
categories night-only and mixed due to data availability.

We placed $\mathrm{N}(0,1)$ priors on each of the $\beta$’s and $\delta$’s,
independent $\mathrm{N}(0,100)$ priors on the off diagonal elements of
$\mathbf{A}$, $\mathrm{Cauchy}^{+}(0,1)$ on each of the range parameters
$\phi_{i}$ and diagonal elements of $\mathbf{A}$, and a flat prior on the cut
points. Ordering of the cut points in the posterior is imposed by the
likelihood of the latent variables. For each land-use filtering model, we ran
each MCMC for 200,000 interations discarding the first 100,000 as burn in.

###  4.1 Model comparison

We compare the land-use filtering model for each collective efficacy component
to a model that assumes a single stationary latent process. The only
difference in the models is the specification of the correlation matrix. As
such, we use the same $\mathrm{N}(0,1)$ priors on each of the $\beta$’s and
$\delta$’s and a flat prior on the cut points. While the land-use filtering
model defines the correlation matrix,
$\boldsymbol{R}(\boldsymbol{\theta}^{*})$, as in Equation 3, the stationary
model defines the correlation as
$\boldsymbol{R}(\bar{\boldsymbol{\theta}})=(1-\kappa)\boldsymbol{R}(\phi)+\kappa\mathbf{I}$,
where $\kappa=1/(1+\tau^{2})$, $\mathbf{I}$ is the identity matrix,
$\boldsymbol{R}(\phi)$ is defined using the exponential correlation function
and $\bar{\boldsymbol{\theta}}=\\{\phi,\tau^{2}\\}$. We place
$\mathrm{Cauchy}^{+}(0,1)$ priors on $\phi$ and $\tau^{2}$. As we are not
interested in using predictions from the stationary model and the covariance
parameters (i.e. only $\phi$ and $\tau^{2}$) converge more quickly, we obtain
10,000 samples and discard the first 2,000 as burn in.

To compare the fitted models, we use the Widely Applicable Information
Criterion (WAIC) (Watanabe, 2010). Information criteria use the log predictive
density and include a penalty for over fitting. In our models, the predictive
density (or likelihood) for the data augmentation model is recovered by
integrating out the latent variable $\mathbf{Z}^{*}$, resulting in an integral
over the multivariate normal distribution:

$p(\mathbf{Y}|\boldsymbol{\theta}^{*},\boldsymbol{\delta},\boldsymbol{\beta},\boldsymbol{\gamma})=\int\dots\int
p(\mathbf{Y},\mathbf{Z}^{*}|\boldsymbol{\theta}^{*},\boldsymbol{\delta},\boldsymbol{\beta},\boldsymbol{\gamma})d\mathbf{Z}^{*}.$

WAIC is defined on the log pointwise posterior predictive density, which
requires the posterior predictive density to be factored by individual
observations. In the land-use filtered model this factorization is obtained by
conditioning on the spatial random effects $\tilde{\mathbf{Z}}$ in the model,

$\displaystyle
p(\mathbf{Y}|\boldsymbol{\theta}^{*},\tilde{\mathbf{Z}},\boldsymbol{\beta},\boldsymbol{\gamma})$
$\displaystyle=\int\dots\int
p(\mathbf{Y},\mathbf{Z}^{*}|\boldsymbol{\theta}^{*},\tilde{\mathbf{Z}},\boldsymbol{\beta},\boldsymbol{\gamma})d\mathbf{Z}^{*}$
$\displaystyle=\int_{C_{1}}p(Z_{1}^{*}|\boldsymbol{\theta}^{*},\tilde{\mathbf{Z}},\boldsymbol{\beta},\boldsymbol{\gamma})dZ_{1}^{*}\dots\int_{C_{n}}p(Z_{n}^{*}|\boldsymbol{\theta}^{*},\tilde{\mathbf{Z}},\boldsymbol{\beta},\boldsymbol{\gamma})dZ_{n}^{*}$
$\displaystyle=\prod_{r=1}^{n}[\Phi(\gamma_{Y_{r}}-\mathbf{x}_{r}^{\prime}\boldsymbol{\beta}-\mathbf{h}_{r}^{\prime}\tilde{\mathbf{Z}}))-\Phi(\gamma_{Y_{r}-1}-\mathbf{x}_{r}^{\prime}\boldsymbol{\beta}-\mathbf{h}_{r}^{\prime}\tilde{\mathbf{Z}})],$

where $\Phi()$ is the standard normal cumulative distribution function and
$\mathbf{h}_{r}$ is the column vector of the $r$th row of $\mathbf{H}$. The
distribution of $Y_{r}$ given $\mathrm{Z}^{*}_{r}$ is simply the indicator
function
$\mathrm{I}(\boldsymbol{\gamma}_{Y_{r}-1}<\mathrm{Z}^{*}_{r}<\boldsymbol{\gamma}_{Y_{r}})$,
which gives the above bounds of integration

$C_{r}=[\gamma_{Y_{r}-1},\gamma_{Y_{r}}],$

where we define $\gamma_{0}=-\infty$ and $\gamma_{K}=\infty$. As a post
processing step of the MCMC output, we obtain samples from the posterior
distribution of the spatial random effects by drawing samples from the full
conditional distribution given in Equation 4 for all observed locations.

###  4.2 Results

We begin with a comparison of model fit using WAIC. The WAIC of the land-use
filtering model and stationary model for each component of collective efficacy
is given in Table 3. To calculate WAIC we drew the spatial random effect from
thinned (every 8th draw after burn in for the stationary models and every
100th draw after burn in for the land-use filtering models) MCMC chains to
give us 1000 draws from each model. Thinning facilitated a quicker computation
of WAIC. For all three components, the WAIC is lower for the land-use
filtering model, indicating a better model fit. The dimension expansion
technique applied in the land-use filtering model allows us to estimate
different spatial processes for each land-use type, allowing for a better fit
of the underlying spatial process for each component of collective efficacy
across land-use types.

Table 3: WAIC scores for the land-use filtering and the stationary model | Defense |  | Trust |  | Observation
---|---|---|---|---|---
Land-use Filtering | 17013.3 |  | 13203.9 |  | 17998.1
Stationary Covariance | 17732.3 |  | 13968.6 |  | 18871.0

Now we discuss how parameter estimates from our land-use filtering and
stationary models compare, noting that these comparisons should be viewed
heuristically given the partial identifiability of the covariance parameters
characterizing the strength of spatial dependence. Tables 4 and 5 give the
posterior mean estimates of the parameters in the land-use filtering and
stationary models respectively. For the stationary models, the single estimate
of $\phi$ is in between estimates of $\phi_{1}$ and $\phi_{2}$ from the land-
use filtering model. The stationary model smooths evenly over the entire study
region and masks variation by land-use category.

Table 4: Posterior mean estimates and corresponding posterior 95% credible intervals of the parameters in the land-use filtered Bayesian ordinal regression model. The numeric subscripts for the $\beta$’s and $\delta$’s correspond to the land-use processes: 1 for residential, 2 for nonresidential, 3 for the “other” category. Note, with the LMC, the subscripts for the $a_{qq^{\prime}}$’s correspond to the linear weights that relate the independent component processes $w_{q}(\mathbf{s})$, each with spatial range $\phi_{q}$, to the multivariate outcome $\boldsymbol{\eta}(\mathbf{s})$. For example, the residential process is defined by $a_{11}$ and $\phi_{1}$, the nonresidential process is defined by $a_{21}$, $a_{22}$, $\phi_{1}$ and $\phi_{2}$. The third component process $w_{3}(\mathbf{s})$ is spatially independent, thus $\phi_{3}$ is fixed at positive infinity. The subscripts “tod” and “dow” stand for time of day and day of week with “n” for night, “m” for mixed, and “e” for weekend. | Defense |  | Trust |  | Observation
---|---|---|---|---|---
|  Posterior mean |  95% Credible interval |  |  Posterior mean |  95% Credible interval |  |  Posterior mean |  95% Credible interval
$\phi_{1}$ | 325.06 | (259, | 408) |  | 405.40 | (333, | 486) |  | 2142.49 | (1215, | 4797)
$\phi_{2}$ | 2.54 | (0.43, | 7.88) |  | 24.29 | (11.5, | 41.7) |  | 92.70 | (58.6, | 136.0)
$a_{11}$ | 1.39 | (1.22, | 1.56) |  | 1.76 | (1.61, | 1.92) |  | 1.60 | (1.48, | 1.73)
$a_{21}$ | 0.32 | (0.22, | 0.41) |  | 0.31 | (0.18, | 0.43) |  | 0.05 | (-0.12, | 0.20)
$a_{22}$ | 0.58 | (0.42, | 0.72) |  | 0.68 | (0.56, | 0.81) |  | 0.34 | (0.28, | 0.39)
$a_{31}$ | 0.64 | (0.02, | 2.39) |  | 0.79 | (0.06, | 1.89) |  | 0.53 | (0.03, | 1.44)
$a_{32}$ | -0.39 | (-3.91, | 3.56) |  | 0.71 | (-0.14, | 1.64) |  | 0.63 | (-0.07, | 1.45)
$a_{33}$ | 2.42 | (0.15, | 6.89) |  | 0.67 | (0.03, | 1.89) |  | 0.64 | (0.03, | 1.76)
$\delta_{2}$ | 0.08 | (-0.62, | 0.77) |  | -0.37 | (-0.74, | 0.02) |  | 0.22 | (0.10, | 0.33)
$\delta_{3}$ | 0.03 | (-0.76, | 0.80) |  | -0.34 | (-0.81, | 0.15) |  | 0.10 | (-0.52, | 0.72)
$\beta_{1,\mathrm{tod=n}}$ | -0.40 | (-0.46, | -0.35) |  | -0.20 | (-0.25, | -0.15) |  | -0.18 | (-0.23, | -0.13)
$\beta_{1,\mathrm{tod=m}}$ | 0.12 | (-0.00 | 0.25) |  | -0.12 | (-0.26, | 0.02) |  | 0.15 | (0.01, | 0.29)
$\beta_{2,\mathrm{tod=n}}$ | -0.02 | (-0.11, | 0.06) |  | -0.07 | (-0.17, | 0.02) |  | -0.12 | (-0.21, | -0.03)
$\beta_{2,\mathrm{tod=m}}$ | 0.01 | (-0.07, | 0.09) |  | -0.04 | (-0.13, | 0.06) |  | 0.01 | (-0.08, | 0.11)
$\beta_{2,\mathrm{dow=e}}$ | -0.06 | (-0.14, | 0.02) |  | -0.02 | (-0.11, | 0.08) |  | -0.23 | (-0.32, | -0.14)
$\beta_{2,\mathrm{dow=m}}$ | -0.04 | (-0.11, | 0.04) |  | -0.01 | (-0.10, | 0.07) |  | -0.11 | (-0.19, | -0.03)
$\beta_{3,\mathrm{tod=n/m}}$ | -0.17 | (-0.75, | 0.42) |  | -0.03 | (-0.73, | 0.70) |  | 0.04 | (-0.59, | 0.66)
$\beta_{3,\mathrm{dow=e/m}}$ | -0.07 | (-0.64, | 0.50) |  | -0.15 | (-0.80, | 0.50) |  | -0.34 | (-0.94, | 0.25)
$\gamma_{1}$ | -1.84 | (-1.94, | -1.75) |  | -1.98 | (-2.07, | -1.88) |  | -1.93 | (-2.00, | -1.87)
$\gamma_{2}$ | -1.25 | (-1.33, | -1.17) |  | -1.33 | (-1.40, | -1.27) |  | -1.26 | (-1.32, | -1.20)
$\gamma_{3}$ | -0.50 | (-0.57, | -0.43) |  | -0.42 | (-0.48, | -0.37) |  | -0.25 | (-0.31, | -0.20)
$\gamma_{4}$ | 0.72 | (0.64, | 0.80) |  | 0.72 | (0.66, | 0.79) |  | 0.65 | (0.60, | 0.70)

The regression coefficients for the two models are comparable. The baseline
category is residential locations during the daytime. For all three
components, we estimate a lower level in the latent process (i.e. higher
probability of a lower rating) for residential locations at night. We are able
to estimate this effect as caregivers were asked to rate their home
neighborhood on all three components of collective efficacy for both the
daytime and the nighttime. For nonresidential and “other” locations, we do not
explicitly have this counter factual information in the data set. For the
land-use mean shift parameters in the land-use filtering model, we estimate an
increase in observation for nonresidential locations compared to residential
locations. Additionally, there is some evidence of an negative effect on trust
for nonresidential locations compared to residential locations. The stationary
covariance model estimates more clear effects (credible intervals that do not
contain zero) of the mean shift parameters which may arise due to over-
smoothing across land-use categories.

Table 5: Estimated posterior means and corresponding posterior 95% credible intervals of the parameters in the stationary Bayesian ordinal regression model. The numeric subscripts for the $\beta$’s and $\delta$’s correspond to the land-use processes: 1 for residential, 2 for nonresidential, 3 for the “other” category. The subscripts “tod” and “dow” stand for time of day and day of week with “n” for night, “m” for mixed, and “e” for weekend. | Defense |  | Trust |  | Observation
---|---|---|---|---|---
|  Posterior mean |  95% Credible interval |  |  Posterior mean |  95% Credible interval |  |  Posterior mean |  95% Credible interval
$\phi$ | 103.21 | (69.8, | 140.0) |  | 120.15 | (86.5, | 155.1) |  | 494.95 | (348.4, | 716.5)
$\tau^{2}$ | 0.33 | (0.26, | 0.42) |  | 0.79 | (0.66, | 0.92) |  | 0.31 | (0.25, | 0.38)
$\delta_{2}$ | 0.09 | (0.01, | 0.16) |  | -0.33 | (-0.41, | -0.24) |  | 0.16 | (0.09, | 0.23)
$\delta_{3}$ | 0.09 | (-0.24, | 0.41) |  | -0.22 | (-0.58, | 0.15) |  | 0.09 | (-0.54, | 0.71)
$\beta_{1,\mathrm{tod=n}}$ | -0.41 | (-0.48, | -0.34) |  | -0.16 | (-0.22, | -0.09) |  | -0.24 | (-0.30, | -0.18)
$\beta_{1,\mathrm{tod=m}}$ | 0.13 | (-0.01, | 0.26) |  | -0.09 | (-0.23, | 0.06) |  | 0.09 | (-0.05, | 0.23)
$\beta_{2,\mathrm{tod=n}}$ | -0.04 | (-0.12, | 0.05) |  | -0.08 | (-0.17, | 0.01) |  | -0.12 | (-0.21, | -0.03)
$\beta_{2,\mathrm{tod=m}}$ | 0.01 | (-0.08, | 0.09) |  | -0.03 | (-0.12, | 0.06) |  | 0.01 | (-0.08, | 0.10)
$\beta_{2,\mathrm{dow=e}}$ | -0.06 | (-0.14, | 0.03) |  | -0.02 | (-0.12, | 0.07) |  | -0.23 | (-0.32, | -0.14)
$\beta_{2,\mathrm{dow=m}}$ | -0.00 | (-0.08, | 0.07) |  | 0.00 | (-0.08, | 0.09) |  | -0.10 | (-0.18, | -0.02)
$\beta_{3,\mathrm{tod=n/m}}$ | -0.10 | (-0.76, | 0.55) |  | -0.13 | (-0.83, | 0.59) |  | 0.02 | (-0.62, | 0.66)
$\beta_{3,\mathrm{dow=e/m}}$ | 0.03 | (-0.57, | 0.64) |  | -0.17 | (-0.80, | 0.45) |  | -0.40 | (-1.00, | 0.20)
$\gamma_{1}$ | -1.97 | (-2.02, | -1.92) |  | -1.77 | (-1.86, | -1.67) |  | -2.00 | (-2.04, | -1.95)
$\gamma_{2}$ | -1.34 | (-1.37, | -1.30) |  | -1.13 | (-1.20, | -1.06) |  | -1.34 | (-1.37, | -1.30)
$\gamma_{3}$ | -0.55 | (-0.57, | -0.53) |  | -0.24 | (-0.31, | -0.16) |  | -0.30 | (-0.33, | -0.29)
$\gamma_{4}$ | 0.75 | (0.72, | 0.78) |  | 0.90 | (0.84, | 0.95) |  | 0.59 | (0.57, | 0.61)

###  4.3 Predictions

Visualization of the predicted levels of each component of collective efficacy
also supports the use of the land-use filtering model. Figure 5 contains the
point-wise posterior mean predictions of the spatial random effect
$\tilde{\mathbf{Z}}$ across a fine grid of points (around 30,000) within the
I-270 belt loop. The land-use category of each prediction location was
assigned by the land-use category of the nearest parcel. We obtained
predictions using Equation 4 for each of the draws of the thinned MCMC chains
used in the WAIC calculations. Plotted is the point-wise posterior mean of the
spatial random effect ($\boldsymbol{\mu}^{[b]}$ in Equation 4) at each
prediction location averaged across the thinned draws. We calculated the
posterior mean of the spatial random effect rather than the mean of posterior
predictive distribution due to memory constraints for obtaining a multivariate
normal sample at a very fine grid of locations. The fine grid of locations
allows us to better visualize within neighborhood variation across land-use
categories. The predicted stationary process on the right is more smooth
compared to the predicted land-use process on the left. Land-use filtering
allows use to estimate distinct processes for each land-use type while also
smoothing over land-use boundaries. There are discontinuities in the
prediction at the boundaries from one land-use type to another, however, these
discontinuities would be greater if we simply fit three independent process
(i.e. the matrix of weights $\mathbf{A}$ is a diagonal matrix) for each land-
use category and then filtered the processes. The predicted land-use filtered
process for the defense and trust components shows the shorter range spatial
dependence within residential areas and the longer range spatial dependence
for non-residential locations. For the observation component, the spatial
range of residential locations is effectively zero, thus the predictions
revert to the mean for residential locations while allowing for greater
spatial smoothing over nonresidential locations.

Figure 5: Point-wise posterior means of the spatial random effects for defense
(top row), trust (middle row), and observation (bottom row). On the left are
the predictions of the collective efficacy components derived from the land-
use filtering model. Predictions from the stationary model are on the right.
The scales are the same for plots on the same row.

As the full Bayesian spatial ordinal regression model is time consuming to fit
we also compared the predictions from the ordinal land-use filtering model to
predictions from the approximate land-use filtering model described in Section
3.5.2. The approximate model treats the rating response directly as a Gaussian
outcome, thus the scales of predictions are different. Additionally, treating
the outcome directly as Gaussian allows one to drop the restrictions on the
covariance function when defining the land-use filtering process. Predictions
from both models are very highly correlated for each component of collective
efficacy: predictions for the defense component had a correlation of 0.93, for
trust a correlation of 0.90, and for observation a correlation of 0.92. In
practice, when uncertainty estimates of the collective efficacy process are
not needed, using the approximate model will give nearly the same results as
the full Bayesian spatial ordinal model with a significant reduction in
computation time.

##  5 Simulation Study

To further understand the uses and limitations of our model we explore model
performance when data are generated from a land-use filtering process as
described in Section 3. In this model generated data setting we seek to
understand the limits of our partially identified model.

As we have a partially identified model, the constraints explained in Section
3.3 do not guarantee that we will be able to recover the values of model
parameters used to generate the data. Indeed, all spatial probit regression
models are partially identifiable: a direct result of binning the latent
spatial process into ordinal or binary categories. While this loss of
information weakens identifiability of the spatial parameters, we can still
evaluate the predictive performance of our model by comparing predictions from
the model to the true latent process of the generative model over a fine grid
of points. This evaluation of model performance aligns with the main goal of
our analysis to understand the fine grain variation of collective efficacy
within census units.

For our simulation study we use the same parameter settings as in the
demonstration of the generative model in section 3.1. Our test set to evaluate
predictive performance is the same fine grid of ten thousand points over the
unit square domain. For each iteration of the simulation study we randomly
generate a training set of one thousand $x$ and $y$ coordinates drawn
uniformly from the unit square. Each training set location is assigned the
land use category of the nearest test set location. Next, we draw a
realization of the latent variable $\mathbf{Z}^{*}$ at the combined set of
test and training locations with mean zero and covariance function defined by
equation 3 with

$\mathbf{A}=\begin{bmatrix}1.8&0.0&0.0\\\ 0.8&1.2&0.0\\\
0.9&1.0&1.25\end{bmatrix},\;\;\boldsymbol{\phi}=\begin{bmatrix}40.0&10.0&2.0\end{bmatrix}^{\prime},\text{
and
}\boldsymbol{\sigma}^{2}=\begin{bmatrix}1.0&1.0&1.0\end{bmatrix}^{\prime}.$

From the latent process realization we apply the cutoffs
$\boldsymbol{\gamma}=\begin{bmatrix}-2.0&-0.5&0.1&1.1\end{bmatrix}^{\prime}$
to obtain the ordinal response. Then we estimate the parameters of the
approximate model as described in section 3.5.2 for both the land-use
filtering and stationary specification (note that for both models
$\boldsymbol{\delta}$ and $\boldsymbol{\beta}$ are set to be zero). Using the
approximate maximum a posteriori estimates we the generate predictions for the
test set conditional on the observed ordinal responses. We use the approximate
model for computational efficiency and because the approximate model
predictions are highly correlated with the full Bayesian model predictions.
Since the predictions conditional on the ordinal data will be on a different
scale than the latent process, to evaluate model performance we compare the
centered and scaled predictions given the ordinal response to the true latent
process values at the test set of locations. Out of 100 replications of the
simulation study (generating data from the land-use filtering model and
fitting both the approximate land-use model and stationary model to the
simulated data at the training locations) the land-use filtering model had a
lower mean (over the test set of locations) absolute error in 76% of the
replications and higher correlation coefficient between the predictions and
true values at the test locations in 81% of the replications. By both the mean
absolute error and correlation coefficient, the land-use filtering model
consistently performs better than the stationary spatial model in predicting
the components of collective efficacy at unobserved locations even though the
model is only partially identifiable.

##  6 Discussion

As demonstrated in the analysis section, our proposed land-use filtering model
for the AHDC collective efficacy data better fits the observed data than the
default approach using a spatial generalized linear model with a single
stationary latent process. This superiority holds for all three components of
collective efficacy. Due to the lack of software for fitting alternative
models with a single nonstationary latent process, as opposed to the default
latent stationary model discussed in Section 4, our model comparison exercise
is somewhat limited in scope. For example, Risser and Calder (2015)’s
covariance-regression model could be used to allow the parameters of the
latent spatial process to vary smoothly with spatially-referenced land-use
covariate information.

In our analysis, we fit our model separately to the ratings of each component
of collective efficacy. While one could consider modeling the ratings of the
three components jointly as a multivariate response, in our opinion, doing so
would limit the applicability of our findings in addressing important
sociological questions. Collective efficacy was originally conceived as a
multidimensional construct comprised of separately operating social processes
(Sampson, Morenoff and Earls, 1999). Modern studies of crime and other
outcomes have looked at the effects of separate components of collective
efficacy (Hipp, 2016; Wickes et al., 2019). Therefore, as an input to
downstream analyses, separate predictions of the component processes will
allow for a richer and more nuanced collection of questions about collective
efficacy to be addressed.

As we had very few reports on locations that fell into the “other” category,
we choose to fix $\phi_{3}=\infty$ the spatial dependence parameter of the
third latent process in the LMC construction. This had the effect of adding
additional independent error to the spatial process for the components of
collective efficacy at the “other” locations. Alternatively, one could
consider dropping the third latent process entirely by fixing $a_{33}=0$, so
that the spatial process of our three land-use categories is an expansion of a
two dimensional latent spatial process. We choose to keep $\phi_{3}$ in the
model to represent our belief that there truly are three underlying, land-use-
specific processes for the components of collective efficacy, and we fix
$\phi_{3}$ due to data availability.

A limitation of our analysis is that we were not able to account for
dependence in the collective efficacy component ratings made by the same
individual. Conceptually, it may be possible to account for differences in how
individuals use the component response scales since individuals report on
several locations, which may be rated by multiple individuals. See Linero,
Bradley and Desai (2018) for one approach to address individual rater effects.
However, since many routine locations will be closer to the caregiver’s home
location, we found that an individual random effect is essentially confounded
with the spatial random effect. Since the latter is of primary interest and
necessary for spatial prediction at unrated locations, we decided to forego an
individual-level random effect in our analyses.

Finally, we note that our prediction of the collective efficacy components
across the study area implicitly assume that the caregivers who rate locations
other than their residence are a representative sample of individuals who
frequent the locations as part of their daily routine. We view this assumption
as foundational to our notion of collective efficacy as a continuously-indexed
spatial process. Alternatively, one could define collective efficacy of a
location from the perspective of an objective observer, who may or may not
spend time at the location. This notion is not compatible with the AHDC Study
design and, more importantly, is arguably not consistent with sociological
theory which emphasizes perceptions of space among users of the space.

In future work, we will use the estimated continuously-indexed spatial
collective efficacy component processes to examine the relationship between
collective efficacy and point-level crime data. We will also explore the
hypothesis that within-neighborhood variation in collective efficacy relates
to crime (Weisburd et al., 2016). Such an analysis is only possible with our
fine-grained maps of collective efficacy across the study area.

In conclusion, while our proposed land-use filtering model is a novel method
for learning about small-scale variability in collective efficacy across an
urban environment, we acknowledge that the approach may be more generally
applicable. Other spatial prediction problems in which the spatial dependence
structure may depend on land use (e.g., pollution mapping) might benefit from
the modeling approach. More broadly, any spatially-dependent, categorical
variable that defines a partition of a well-defined study area could be used
instead of land use in the approach.

This study was supported in part by the National Institute on Drug Abuse
(Christopher R. Browning; R01DA032371); the Eunice Kennedy Shriver National
Institute on Child Health and Human Development (Catherine A. Calder;
R01HD088545; John Casterline, The Ohio State University Institute for
Population Research, P2CHD058484; Elizabeth Gershoff, The University of Texas
at Austin Population Research Center, P2CHD-042849); and the W. T. Grant
Foundation.

Example dataset and code AHDC data will be deposited to Inter-university
Consortium for Political and Social Research (ICPSR) in publicly available and
restricted access forms (expected in Summer 2023). To ensure participant
privacy and maintenance of data confidentiality, the ADHC caregiver location
reports needed to reproduce the analyses in this paper will only be available
in the restricted access version of the data. Qualified researchers will be
able to submit an application to ICPSR for access to the data. In our
supplementary material, we instead provide a synthetic data set generated from
the posterior predictive distribution of the trust component at randomly
selected locations out of the grid of points used for our predictive maps.
Accompanying the example data set is code to reproduce the analyses done in
the paper.

## References

  * (2021) [author] , Franklin County AuditorF. C. A. (2021). Franklin county GIS parcel shape files. 2014 Data retrieved from https://franklincountyauditor.com/ftp.
  * Albert and Chib (1993) [author] Albert, James H.J. H. and Chib, SiddharthaS. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association 88 669–679.
  * Albert and Chib (1997) [author] Albert, JimJ. and Chib, SiddharthaS. (1997). Bayesian methods for cumulative, sequential and two-step ordinal data regression models.
  * Apanasovich, Genton and Sun (2012) [author] Apanasovich, Tatiyana VT. V., Genton, Marc GM. G. and Sun, YingY. (2012). A valid Matérn class of cross-covariance functions for multivariate random fields with any number of components. Journal of the American Statistical Association 107 180-193.
  * Bandura (1982) [author] Bandura, AlbertA. (1982). Self-efficacy mechanism in human agency. The American psychologist 37 122-147.
  * Bandura (1986) [author] Bandura, AlbertA. (1986). Social foundations of thought and action : a social cognitive theory / Albert Bandura. Prentice-Hall series in social learning theory. Prentice-Hall, Englewood Cliffs, N.J.
  * Banerjee, Carlin and Gelfand (2014) [author] Banerjee, SudiptoS., Carlin, Bradley P.B. P. and Gelfand, Alan E.A. E. (2014). Hierarchical modeling and analysis for spatial data, Second ed. Chapman and Hall/CRC, New York. 10.1201/b17115
  * Berrett and Calder (2012) [author] Berrett, CandaceC. and Calder, Catherine A.C. A. (2012). Data augmentation strategies for the Bayesian spatial probit regression model. Computational Statistics & Data Analysis 56 478-490.
  * Berrett and Calder (2016) [author] Berrett, CandaceC. and Calder, Catherine A.C. A. (2016). Bayesian spatial binary classification. Spatial Statistics 16 72-102. https://doi.org/10.1016/j.spasta.2016.01.004
  * Besançon et al. (2021) [author] Besançon, MathieuM., Papamarkou, TheodoreT., Anthoff, DavidD., Arslan, AlexA., Byrne, SimonS., Lin, DahuaD. and Pearson, JohnJ. (2021). Distributions.jl: Definition and modeling of probability distributions in the JuliaStats ecosystem. Journal of Statistical Software 98 1–30. 10.18637/jss.v098.i16
  * Bezanson et al. (2017) [author] Bezanson, JeffJ., Edelman, AlanA., Karpinski, StefanS. and Shah, Viral BV. B. (2017). Julia: A fresh approach to numerical computing. SIAM review 59 65–98.
  * Bornn, Shaddick and Zidek (2012) [author] Bornn, LukeL., Shaddick, GavinG. and Zidek, James VJ. (2012). Modeling non-stationary processes through dimension expansion. Journal of the American Statistical Association 107 281–289. 10.1080/01621459.2011.646919
  * Browning (2002) [author] Browning, Christopher RC. R. (2002). The span of collective efficacy: Extending social disorganization theory to partner violence. Journal of marriage and family 64 833-850.
  * Calder (2008) [author] Calder, Catherine A.C. A. (2008). A dynamic process convolution approach to modeling ambient particulate matter concentrations. Environmetrics 19\. 10.1002/env.852
  * Cohen et al. (2006) [author] Cohen, Deborah AD. A., Finch, Brian KB. K., Bower, AimeeA. and Sastry, NarayanN. (2006). Collective efficacy and obesity: The potential influence of social factors on health. Social science & medicine (1982) 62 769-778.
  * Coleman (1988) [author] Coleman, James SJ. S. (1988). Social capital in the creation of human capital. American journal of sociology 94 S95–S120.
  * Corcoran et al. (2018) [author] Corcoran, JonathanJ., Zahnow, ReneeR., Wickes, RebeccaR. and Hipp, JohnJ. (2018). Neighbourhood land use features, collective efficacy and local civic actions. Urban studies (Edinburgh, Scotland) 55 2372-2390.
  * Cowles (1996) [author] Cowles, Mary KathrynM. K. (1996). Accelerating Monte Carlo Markov chain convergence for cumulative-link generalized linear models. Statistics and Computing 6 101–111.
  * De Oliveira (1997) [author] De Oliveira, VictorV. (1997). Prediction in some classes of non-Gaussian random fields, PhD thesis, University of Maryland, College Park.
  * De Oliveira (2000) [author] De Oliveira, VictorV. (2000). Bayesian prediction of clipped Gaussian random fields. Computational Statistics & Data Analysis 34 299-314.
  * Frigo and Johnson (2005) [author] Frigo, MatteoM. and Johnson, Steven G.S. G. (2005). The design and implementation of FFTW3. Proceedings of the IEEE 93 216–231. Special issue on “Program Generation, Optimization, and Platform Adaptation”. 10.1109/JPROC.2004.840301
  * Gelfand et al. (2004) [author] Gelfand, Alan EA. E., Schmidt, Alexandra MA. M., Banerjee, SudiptoS. and Sirmans, C FC. F. (2004). Nonstationary multivariate process modeling through spatially varying coregionalization. Test (Madrid, Spain) 13 263-312.
  * Gneiting, Kleiber and Schlather (2010) [author] Gneiting, TilmannT., Kleiber, WilliamW. and Schlather, MartinM. (2010). Matérn cross-covariance functions for multivariate random fields. Journal of the American Statistical Association 105 1167-1177.
  * Goulard and Voltz (1992) [author] Goulard, MM. and Voltz, MM. (1992). Linear coregionalization model: Tools for estimation and choice of cross-variogram matrix. Mathematical geology 24 269-286.
  * Grzebyk and Wackernagel (1994) [author] Grzebyk, MichelM. and Wackernagel, HansH. (1994). Multivariate analysis and spatial/temporal scales: Real and complex models. In Proceedings of the XVIIth International Biometrics Conference 1 19–33. Citeseer.
  * Gustafson (2015) [author] Gustafson, PaulP. (2015). Bayesian Inference for Partially Identified Models: Exploring the Limits of Limited Data, 1st ed. Chapman and Hall/CRC. https://doi.org/10.1201/b18308
  * Higgs and Hoeting (2010) [author] Higgs, Megan DaileyM. D. and Hoeting, Jennifer A.J. A. (2010). A clipped latent variable model for spatially correlated ordered categorical data. Comput. Stat. Data Anal. 54 1999–2011. 10.1016/j.csda.2010.02.024
  * Hipp (2016) [author] Hipp, John R.J. R. (2016). Collective efficacy: How is it conceptualized, how is it measured, and does it really matter for understanding perceived neighborhood crime and disorder? Journal of Criminal Justice 46 32-44.
  * Ingebrigtsen, Lindgren and Steinsland (2014) [author] Ingebrigtsen, RikkeR., Lindgren, FinnF. and Steinsland, IngelinI. (2014). Spatial models with explanatory variables in the dependence structure. Spatial Statistics 8 20-38.
  * Jacobs (1961) [author] Jacobs, JaneJ. (1961). The death and life of great American cities. Random House, New York.
  * Li, Ding and Mealli (2022) [author] Li, FanF., Ding, PengP. and Mealli, FabriziaF. (2022). Bayesian Causal Inference: A Critical Review.
  * Lin et al. (2019) [author] Lin, DahuaD., White, John MylesJ. M., Byrne, SimonS., Bates, DouglasD., Noack, AndreasA., Pearson, JohnJ., Arslan, AlexA., Squire, KevinK., Anthoff, DavidD., Papamarkou, TheodoreT., Besançon, MathieuM., Drugowitsch, JanJ., Schauer, MoritzM. and contributors (2019). JuliaStats/Distributions.jl: a Julia package for probability distributions and associated functions. 10.5281/zenodo.2647458
  * Linero, Bradley and Desai (2018) [author] Linero, Antonio R.A. R., Bradley, Jonathan R.J. R. and Desai, ApurvaA. (2018). Multi-rubric models for ordinal spatial data with application to online ratings data. The Annals of Applied Statistics 12 2054 – 2074.
  * Matsueda and Drakulich (2016) [author] Matsueda, Ross LR. L. and Drakulich, Kevin MK. M. (2016). Measuring collective efficacy: A multilevel measurement model for nested data. Sociological methods & research 45 191-230.
  * Mogensen and Riseth (2018) [author] Mogensen, Patrick KofodP. K. and Riseth, Asbjørn NilsenA. N. (2018). Optim: A mathematical optimization package for Julia. Journal of Open Source Software 3 615\. 10.21105/joss.00615
  * Molnar et al. (2015) [author] Molnar, Beth EB. E., Goerge, Robert MR. M., Gilsanz, PaolaP., Hill, AndreaA., Subramanian, S. VS. V., Holton, John KJ. K., Duncan, Dustin TD. T., Beatriz, Elizabeth DE. D. and Beardslee, William RW. R. (2015). Neighborhood-level social processes and substantiated cases of child maltreatment. Child abuse & neglect 51 41-53.
  * Neto, Schmidt and Guttorp (2014) [author] Neto, Joaquim Henriques ViannaJ. H. V., Schmidt, Alexandra M.A. M. and Guttorp, PeterP. (2014). Accounting for spatially varying directional effects in spatial covariance structures. Journal of the Royal Statistical Society. Series C (Applied Statistics) 63 103–122.
  * Pebesma (2018) [author] Pebesma, EdzerE. (2018). Simple features for R: standardized support for spatial vector data. The R Journal 10 439–446. 10.32614/RJ-2018-009
  * Rackauckas and Nie (2017) [author] Rackauckas, ChristopherC. and Nie, QingQ. (2017). DifferentialEquations.jl – A performant and feature-rich ecosystem for solving differential equations in Julia. The Journal of Open Research Software 5\. Exported from https://app.dimensions.ai on 2019/05/05. 10.5334/jors.151
  * Raudenbush and Sampson (1999) [author] Raudenbush, Stephen WS. W. and Sampson, Robert JR. J. (1999). Ecometrics: Toward a science of assessing ecological settings, with application to the systematic social observation of neighborhoods. Sociological methodology 29 1-41.
  * Reich et al. (2011) [author] Reich, Brian JB. J., Eidsvik, JoJ., Guindani, MicheleM., Nail, Amy JA. J. and Schmidt, Alexandra MA. M. (2011). A class of covariate-dependent spatiotemporal covariance functions. The annals of applied statistics 5 2265–2687.
  * Risser and Calder (2015) [author] Risser, Mark D.M. D. and Calder, Catherine A.C. A. (2015). Regression‐based covariance functions for nonstationary spatial modeling. Environmetrics. 26\.
  * Risser et al. (2019) [author] Risser, Mark D.M. D., Calder, Catherine A.C. A., Berrocal, Veronica J.V. J. and Berrett, CandaceC. (2019). Nonstationary spatial prediction of soil organic carbon: Implications for stock assessment decision making. The Annals of Applied Statistics.
  * Sampson, Morenoff and Earls (1999) [author] Sampson, Robert J.R. J., Morenoff, Jeffrey D.J. D. and Earls, FeltonF. (1999). Beyond Social Capital: Spatial Dynamics of Collective Efficacy for Children. American Sociological Review 64 633–660.
  * Sampson, Raudenbush and Earls (1997) [author] Sampson, Robert J.R. J., Raudenbush, Stephen W.S. W. and Earls, FeltonF. (1997). Neighborhoods and violent crime: A multilevel study of collective efficacy. Science 277 918-924. 10.1126/science.277.5328.918
  * Schliep and Hoeting (2015) [author] Schliep, Erin ME. M. and Hoeting, Jennifer AJ. A. (2015). Data augmentation and parameter expansion for independent or spatially correlated ordinal data. Computational statistics & data analysis 90 1-14.
  * Schmidt and Gelfand (2003) [author] Schmidt, Alexandra M.A. M. and Gelfand, Alan E.A. E. (2003). A Bayesian coregionalization approach for multivariate pollutant data. Journal of Geophysical Research: Atmospheres 108\.
  * Schmidt, Guttorp and O’Hagan (2011) [author] Schmidt, Alexandra M.A. M., Guttorp, PeterP. and O’Hagan, AnthonyA. (2011). Considering covariates in the covariance structure of spatial processes. Environmetrics. 22\.
  * Shaw and McKay (1942) [author] Shaw, Clifford RobeC. R. and McKay, Henry DonaldH. D. (1942). Juvenile delinquency and urban areas. University of Chicago Press, Chicago, Ill.
  * R Core Team (2022) [author] R Core Team (2022). R: A language and environment for statistical computing R Foundation for Statistical Computing, Vienna, Austria.
  * Wackernagel (2003) [author] Wackernagel, Hans.H. (2003). Multivariate geostatistics : an introduction with applications / Hans Wackernagel., 3rd, completely rev. ed. ed. Springer, Berlin ;.
  * Watanabe (2010) [author] Watanabe, SumioS. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular Learning Theory. Journal of Machine Learning Research 11 3571-3594.
  * Weisburd et al. (2016) [author] Weisburd, DavidD., Eck, John E.J. E., Braga, Anthony A.A. A., Telep, Cody W.C. W., Cave, BreanneB., Bowers, KateK., Bruinsma, GerbenG., Gill, CharlotteC., Groff, Elizabeth R.E. R., Hibdon, JulieJ. and et al. (2016). Place matters: Criminology for the twenty-first century. Cambridge University Press. 10.1017/CBO9781139342087
  * Wickes et al. (2019) [author] Wickes, RebeccaR., Zahnow, ReneeR., Corcoran, JonathanJ. and Hipp, John RJ. R. (2019). Neighbourhood social conduits and resident social cohesion. Urban studies (Edinburgh, Scotland) 56 226-248.
  * Wickham et al. (2019) [author] Wickham, HadleyH., Averick, MaraM., Bryan, JenniferJ., Chang, WinstonW., McGowan, Lucy D’AgostinoL. D., François, RomainR., Grolemund, GarrettG., Hayes, AlexA., Henry, LionelL., Hester, JimJ., Kuhn, MaxM., Pedersen, Thomas LinT. L., Miller, EvanE., Bache, Stephan MiltonS. M., Müller, KirillK., Ooms, JeroenJ., Robinson, DavidD., Seidel, Dana PaigeD. P., Spinu, VitalieV., Takahashi, KohskeK., Vaughan, DavisD., Wilke, ClausC., Woo, KaraK. and Yutani, HiroakiH. (2019). Welcome to the tidyverse. Journal of Open Source Software 4 1686\. 10.21105/joss.01686
  * Yan et al. (2007) [author] Yan, JunJ., Cowles, Mary KathrynM. K., Wang, ShaowenS. and Armstrong, Marc P.M. P. (2007). Parallelizing MCMC for Bayesian spatiotemporal geostatistical models. Statistics and Computing 17 323–335.