# Local Spatial Autocorrelation (2)

## Introduction

In this chapter, we will explore some more advanced aspects of local spatial autocorrelation statistics, focusing on bivariate and multivariate extensions. We will first cover the local counterparts of the bivariate Moran’s I, the differential Moran’s I and the Empirical Bayes (EB) standardized Moran’s I. Next, we introduct a bivariate and multivariate extension of the local Join Count statistic, for cases where the two variable can coincide in the same location (so-called, co-location), and (less common) cases where they cannot. Finally, we cover the multivariate extension of the Local Geary statistic.

### Objectives

• Identify clusters and outliers with the bivariate extensions of the local Moran statistic

• Correct the local Moran statistic for variance instability in rates

• Assess clusters of binary variables with and without co-location with a local Join Count statistic

• Identify multivariate clusters by means of the local Geary statistic

• Interpret interaction effects

• Interpret significance

#### GeoDa functions covered

• Space > Bivariate Local Moran’s I
• Space > Differential Local Moran’s I
• Space > Local Moran’s I with EB rate
• Space > Bivariate Local Join Count
• Space > Multivariate Local Join Count
• Space > Multivariate Local Geary

## Extensions of Local Moran’s I

### Getting Started

To illustrate the extensions of local Moran’s I, we will follow the examples we used for their global counterparts. We will use the familiar U.S. county homicide data set natregimes with queen contiguity weights in natregimes_q. By means of the Time Editor, we first group the homicide rates hr60 through hr90 into the variable hr, the homicide counts hc60 through hc90 into the variable hc, and the population counts po60 through po90 into po.

At this point, we are ready to proceed.

### Bivariate Local Moran’s I

#### Concept

The treatment of the bivariate local Moran’s I closely follows that of its global counterpart (see also Anselin, Syabri, and Smirnov 2002). In essence, it captures the relationship between the value for one variable at location $$i$$, $$x_i$$, and the average of the neighboring values for another variable, i.e., its spatial lag $$\sum_j w_{ij} y_j$$. Apart from a constant scaling factor (that can be ignored), the statistic is the product of $$x_i$$ with the spatial lag of $$y_i$$ (i.e., $$\sum_j w_{ij}y_j$$), with both variables standardized, such that their means are zero and variances equal one:

$I_{B,i} = c x_i \sum_j w_{ij} y_j,$ where $$w_{ij}$$ are the elements of the spatial weights matrix.

As for its global counterpart, this statistic needs to be interpreted with caution, since it ignores in-situ correlation between the two variables (see the discussion of global bivariate spatial autocorrelation for details).

A special case of the bivariate local Moran statistic is comparing the same variable at two points in time. The most meaningful application is where one variable is for time period $$t$$, say $$z_t$$, and the other variable is for the neighbors in the previous time period, say $$z_{t-1}$$. This is illustrated in the example below.

#### Implementation

The bivariate local Moran functionality is invoked from the Cluster Maps toolbar icon, or from the menu as Space > Bivariate Local Moran’s I, as shown in Figure 1. Figure 1: Bivariate Local Moran from the toolbar

Next, we proceed in the usual fashion by selecting a variable for $$x$$, for example hr(90), and a variable for $$y$$, say hr(80), in the variable selection dialog, shown in Figure 2. Since the variables are time enabled, the proper period is selected by means of the Time drop down list. Also note that the spatial weights matrix has been set to natregimes_q. Figure 2: Bivariate Local Moran variable selection

The next prompt gives the familiar choice of a Significance Map, Cluster Map, and Moran Scatter Plot, with just the Cluster Map selected as the default. With the default setting, the result is based on 999 permutations with a p-value of 0.05, as shown in Figure 3. Figure 3: Bivariate Local Moran cluster map

All the options for randomization, significance filter, saving the results, etc., are the same as for the univariate local Moran and will not be further discussed here.

#### Interpretation

As mentioned, the interpretation of the bivariate local Moran cluster map warrants some caution, since it does not control for the correlation between the two variables at each location (i.e., the correlation between $$x_i$$ and $$y_i$$). To illustrate this issue, consider the space-time case, where one can interpret the results of the high-high and low-low clusters in Figure 3 as locations where high/low values at time $$t$$ (in our example, homicide rates in 1990) were surrounded by high/low values at time $$t -1$$ (homicide rates in 1980) more so than would be the case randomly. This could either be attributed to the surrounding locations in the past affecting the central location at the current time (a form of a diffusion effect), but could also be due to a strong inertia, in which case there would be no diffusion. More precisely, if the homicide rates in all locations are highly correlated over time, then if the surrounding values in $$t - 1$$ were correlated with the value at $$i$$ in $$t - 1$$, then they would also tend to be correlated with the value at $$i$$ in $$t$$ (through the correlation between $$y_{j,t-1}$$ and $$y_{j,t}$$). Hence, while these findings could be compatible with diffusion, this is not necessarily the case. The same complication affects the interpretation of the bivariate local Moran coefficient between two variables at the same point in time.

To illustrate the problems with interpreting bivariate local cluster maps, we take the case of spatial outliers. In addition to the bivariate results, we also consider the univariate local Moran cluster maps for the variable in each time period, i.e., in $$t - 1$$ and in $$t$$, i.e., hr(80) and hr(90), as shown in Figures 4 and 5. Figure 4: Local Moran cluster map – hr80 Figure 5: Local Moran cluster map – hr90

The highlighted county is Inyo county (CA), which is identified as a low-high spatial outlier in Figure 3. In 1980, this county was part of a high-high cluster, whereas in 1990, it was not significant. In this instance, we can interpret the spatial outlier in the space-time cluster map as indicating a (relative) decrease in the homicide rate for that county in 1990 relative to 1980, when it was surrounded by high homicide rate counties. However, what we cannot conclude is that the surrounding counties still have high homicide rates in 1990. In fact, the pattern is insufficient to be considered significant.

### Differential Local Moran’s I

#### Concept

The differential local Moran statistic is the local counterpart to the differential Moran scatter plot, discussed in an earlier chapter. Instead of using the observations on a variable at two different time periods, this statistic is based on the change over time, i.e., the difference between $$y_t$$ and $$y_{t-1}$$. Note that this is the actual difference and not the absolute difference, so that a positive change will be viewed as high, and a negative change as low. The differences are used in standardized form, i.e., they are not the differences between the standardized variable at two points in time, but the standardized differences between the original values for the variable.

The formal expression for this statistic follows the same logic as before, and consists of the cross product of the difference between $$y_t$$ and $$y_{t-1}$$ at $$i$$ with the associated spatial lag:

$I_{D,i} = c (y_{i,t} - y_{i,t-1}) \sum_j w_{ij} (y_{j,t} - y_{j,t-1})$ The scaling constant $$c$$ can be ignored. Inference is based on conditional permutation. All the usual caveats hold about multiple comparisons, etc.

#### Implementation

The differential local Moran functionality is invoked from the Cluster Maps toolbar icon, or from the menu as Space > Differential Local Moran’s I, as shown in Figure 6. Figure 6: Differential Local Moran from the toolbar

Next follows a variable selection dialog. This takes on a slightly different form (but the same as for the differential Moran scatter plot), as illustrated in Figure 7. First, the variable of interest is selected (here, hr), and then the two time periods are chosen (here, 90 and 80). Note that the system is agnostic about the actual time periods, so that any combination can be selected. The statistic is computed for the difference between the time period specified as the first item and that given as the second item. In our example, the spatial weights are natregime_q. Figure 7: Differential Local Moran variable selection

As before, a choice is presented between a Significance Map, Cluster Map, and Moran Scatter Plot, with just the Cluster Map selected as the default. With the default setting, the result is based on 999 permutations with a p-value of 0.05, as shown in Figure 8. Figure 8: Differential Local Moran cluster map

All the options for randomization, significance filter, saving the results, etc., are the same as for the univariate local Moran and will not be further discussed here.

#### Interpretation

The first aspect of the results is a much smaller number of significant locations compared to the univariate and bivariate cluster maps. These are locations where the change in the variable over time is matched by similar/dissimilar changes in the surrounding locations. It is important to keep in mind that the focus is on change, and there is no direct connection to whether this changes is from high or low values.

Two situations can be distinguished, depending on whether the change variable takes on both positive and negative values, or when all the changes are of the same sign (i.e., all observations either increase or decrease over time).

When both positive and negative change values are present, the high-high locations will tend to be locations with a large increase (positive change), surrounded by locations with similar large increases. The low-low locations will be observations with a large decrease (negative change), surrounded by locations with similar large decreases. Spatial outliers will be locations where an increase is surrounded by a decrease and vice versa.

When all changes are of the same sign, the interpretation of high-high and low-low depends on the sign. Due to the standardization, large positive values will be considered high (above the mean), whereas large negative values will be labeled low (below the mean). This should be kept in mind when interpreting the results.

#### Saving the results

Similar to the functionality for the differential Moran scatter plot, the Save Results option includes an item to save the actual change variable (in raw form, not in standardized form). In the dialog, this corresponds to the Diff Values item, as shown in Figure 9. The other options are the same as for all local spatial autocorrelation statistics, i.e., the value of the statistic, cluster type, and p-value. Figure 9: Saving the results for differential local Moran

### Local Moran with EB Rate

#### Concept

The last of the extensions of the local Moran’s I pertains to the special case where the variable of interest is a rate or proportion. As discussed for the Moran scatter plot, the resulting variance instability can cause problems for the Moran statistic. The EB standardization suggested by Assunção and Reis (1999) for the global case can be extended to the local statistic in a straightforward manner. The statistic has the usual form, but is computed for the standardized rates, $$z$$.

$I_{EB,i} = c z_i \sum_j w_{ij} z_j,$

The standardization of the raw rate $$r_i$$ is the same as before, and is repeated here for completeness (for a more detailed discussion, see the relevant chapter):

$z_i = \frac{r_i - \beta}{\sqrt{\alpha + (\beta/P_i)}}$ with $$\beta$$ as an estimate of the mean and the denominator as an estimate of the standard error.2

All inference and interpretation is the same as for the univariate case and is not further pursued here.

#### Implementation

The local Moran functionality for standardized rates is invoked from the Cluster Maps toolbar icon, or from the menu as Space > Local Moran’s I with EB Rate, as shown in Figure 10. Figure 10: Local Moran for EB rates from the toolbar

Since the rate standardization is computed as part of the operation, the variable selection interface is similar to that used for rate maps. With the time enabled variables, we need to make sure the Time is synchronized so that both the numerator (the events) and the denominator (the population at risk) pertain to the same period (here, 90). In our example, shown in Figure 11, we take the Event Variable as homicide counts, hc(90), and the Base Variable as population, po(90). As before, we use the queen contiguity, natregimes_q. Figure 11: EB rates local Moran variable selection

Again, we can choose between a Significance Map, a Cluster Map and the Moran Scatter Plot options. With the default settings (cluster map, 999 permutations and p = 0.05), the resulting map is as in Figure 12. Figure 12: Local Moran cluster map for EB rates - hr90

Compared to the cluster map for the raw rates in Figure 5, we observe major differences in the low-low clusters in the upper midwest. Typically, the greater the variation among values for the base variable (population at risk), especially when the latter is small, the more the two maps will tend to differ. However, when the base population is more or less equal by design (e.g., for census tracts in the U.S.), there is little gain from using the EB rates.

#### Saving the results

Here again, the Save Results option includes an item to save the actual EB rate (this is identical to the EB Rate standardization that can be computed with the Table Calculator option). In the dialog, this corresponds to the EB Rates item, as shown in Figure 13. The other options are the same as for all local spatial autocorrelation statistics, i.e., the value of the statistic, cluster type, and p-value. Figure 13: Saving the results for EB rates

## Bivariate and Multivariate Local Join Count Statistics

In addition to introducing a univariate local Join Count statistic, Anselin and Li (2019) also include extensions to a bivariate and multivariate setting. In each instance, the variables under consideration, say $$x$$ and $$z$$, only take on a value of 0 and 1.

We distinguish between two cases, referred to as co-location and no co-location. In the first, it is possible for both variables to take the same value at a location. In orther words, at location $$i$$, it is possible to have $$x_i = z_i = 1$$. For example, this would be the case when the variables pertain to two different types of events that can occur in all locations, such as the presence of several non-exclusive characteristics (e.g., a city block that is both low-density and commercial). In the second instance, whenever $$x_i = 1$$, then $$z_i = 0$$ and vice versa. In other words, the two events cannot happen in the same location. This would be case where the classification of a location is exclusive, i.e., if a location is of one type, it cannot be of another type (e.g., a zoning classification).

In GeoDa, the co-location case is treated under Multivariate Local Join Count (which includes the bivariate setup as a special case), whereas the no co-location case is implemented under Bivariate Local Join Count.

While the natural application of Join Count statistics is to deal with binary variables, they can also be extended to a notion of quantile spatial autocorrelation, which is further elaborated upon below.

### Bivariate - No Co-Location

#### Preliminaries

To illustrate the no co-location case, we will use a data set with population figures for the 77 Community areas in Chicago (commpop).3 This data set contains two binary variables, one with value 1 for areas with positive population change between 2000 and 2010 (popplus, with 17 observations = 1), and the other with value 1 for areas with negative population change (popneg, with 60 observations = 1). In other words, by construction the two variables are the complement of each other so that they cannot both have the same value at the same location.

We illustrate popplus in the unique values map shown in Figure 14.4 Figure 14: Positive population change - Chicago Community Areas

We are interested in instances of negative spatial autocorrelation, where locations with popneg = 1 are surrounded by more locations with popplus = 1 than would be expected under randomness. Note that it is important to formulate the research question in this fashion. Since the majority of locations has popneg = 1 (60/77, or 0.78), the phenomenon of a popplus = 1 to be surrounded by locations with popneg = 1 is not rare. In such a situation, the significance test is not meaningful.5 This will generally be the case when a variable has more than 50% of the observations equal to 1. Since only 0.22 of the observations have popplus = 1, it is of interest to see the extent to which these cluster among themselves (e.g., by means of a univariate local Join Count test), or, group around locations with negative population change. Those instances of popneg = 1 could then be characterized as spatial outliers, i.e., areas with population decline surrounded by neighbors with population increase.

This example is primarily to illustrate the functionality of the bivariate test, and is less of substantive interest. For a more meaningful example, see the empirical illustration in Anselin and Li (2019).

Before we can proceed, we also need a spatial weights matrix. We use queen contiguity in the file commarea_q.gal.6 Note that this weights matrix will be used in unstandardized (binary) form.

#### Concept

In its general form, the expression for a bivariate local join count statistic is $BJC_i = x_i (1 - z_i) \sum_j w_{ij} z_j (1 - x_j),$ with $$w_{ij}$$ as unstandardized (binary) spatial weights.

The roles of $$x$$ and $$z$$ may be reversed, but the statistic is not symmetric, so that the results will tend to differ whether $$x$$ or $$z$$ is the focus. In addition, it is important to keep in mind that the statistic will only be meaningful when the proportion of the second variable (for the neighbors) in the population is small.

Since the condition of no co-location ensures that $$1 - z_i = 0$$ whenever $$x_i = 1$$, and vice versa, the statistic simplifies to: $BJC_i = x_i \sum_j w_{ij} z_j.$

A pseudo p-value can be obtained from a one-sided conditional permutation test. This is implemented by carrying out a series of $$k_i$$ draws (with $$k_i$$ as the number of neighbors) for each location $$i$$ where $$x_i = 1$$ (and thus $$z_i = 0$$). The draws are without replacement from N − 1 data tuples $$(x_j, z_j)$$ of which $$Q$$ observations have $$z = 1$$ (since $$z_i = 0$$) and $$P − 1$$ observations have $$x = 1$$. In practice, we only need to draw the $$z_j$$, since the matching $$x_j$$ are zero by construction. The number of times the resulting local join count statistic equals or exceeds the observed value yields a pseudo p-value, in the usual way.

#### Implementation

The bivariate local Join Count functionality for the no co-location case is invoked from the Cluster Maps toolbar icon, or from the menu as Space > Bivariate Local Join Count, as shown in Figure 15. Figure 15: No Co-Location bivariate local Join Count from the toolbar

Next follows the selection of the two binary variables. The first one ($$x$$) is selected from the left-hand column in the dialog, the second ($$z$$) from the right-hand column. In our example, these are popneg and popplus, as shown in Figure 16. Figure 16: Bivariate local Join Count variable selection

As is the case for the univariate local Join Count statistic, only one map is created, i.e., a significance map that highlights the significant locations. For the default settings with 999 permutation and p = 0.05, the result is as shown in Figure 17, with two significant locations. Figure 17: Bivariate local Join Count significance map

#### Interpretation

As mentioned earlier, in the current context, the result can be interpreted as evidence of negative spatial autocorrelation, i.e., evidence of spatial outliers. We can illustrate this with a unique values map for popneg. Invoking the map option (right click) Connectivity > Show Selection and Neighbors will outline the neighbors of a selected observation. With the selection of the most significant location (the Garfield Ridge community area), its neighbors are highlighted, as shown in Figure 18. Figure 18: Spatial outlier

In our example, the significant area (with negative population change) is surrounded by one other neighbor with negative population change, but four neighbors with positive population change. Having four of the five neighbors with the opposite sign is sufficient to identify this location as a spatial outlier.

#### Saving the results

As is the case for the other local spatial autocorrelation coefficients, the results for each location can be saved to the data table by means of the Save Results option. Three items are included: the join count statistic (i.e., the number of neighbors with value $$z_j = 1$$), the total number of neighbors, and the pseudo p-value, as shown in Figure 19. Figure 19: Bivariate Local Join Count Save Results

In our example, this gives the entries for the two significant locations as in Figure 20. Figure 20: Bivariate Local Join Count results in Table

This confirms that for Garfield Ridge, 4 of the 5 neighbors were of the opposite sign, yielding a pseudo p-value of 0.008.

### Bivariate and Multivariate - Co-Location

#### Preliminaries

We will confine our illustration of this statistic to the bivariate case, since the multivariate setting follows as a straightforward extension. Again, we will use a somewhat contrived example, converting continuous variables to categories. We will take observations in the top quintile for two variables in the Guerry sample data set and turn them into a binary variable. Specifically, we will select the observations in the top quintile for the variables Infants (children born out of wedlock)7 and Donatns (donations). We employ the usual procedure by selecting the category in the quintile map (click on the corresponding box in the legend) and saving the selection. The spatial pattern of the resulting binary variables is illustrated in Figures 21 and 22. Figure 21: Infants - fifth quintile Figure 22: Donations - fifth quintile

We are interested in detecting instances where the co-location of a top quintile for both variables coincides with such a co-location for the neighbors, more so than expected under spatial randomness. In other words, we assess whether a location with high values for both variables is also surrounded by neighboring locations with high values for both.

We can gain an initial insight into the potential for such clusters from a simple co-location map for the two variables, shown in Figure 23. In all, there are only 5 potential significant cluster locations. Figure 23: Co-Location map - top quintile Infants and Donations

As the spatial weights matrix, we use queen contiguity, guerry_85_q. We are now ready to proceed.

#### Concept

A bivariate co-location cluster requires that $$x_i = z_i = 1$$ coincides with neighbors for which also $$x_j = z_j = 1$$. For the bivariate case, the corresponding local Join Count statistic takes the form: $CLC_i = x_iz_i \sum_j w_{ij}x_jz_j,$ with $$w_{ij}$$ as unstandardized (binary) spatial weights.

As before, a conditional permutation approach can be constructed for those locations with $$x_i = z_i = 1$$. We draw $$k_i$$ pairs of observations $$(x_j, z_j)$$ from the set of $$N−1$$ (this contains $$P−1$$ observations with $$x_j = 1$$ and $$Q−1$$ observations with $$z_j = 1$$). In a one-sided test, we again count the number of times the statistic equals or exceeds the observed join count value at $$i$$.

The extension to more than two variables is mathematically straightforward. We consider $$m$$ variables at each location $$i$$, i.e., $$x_{hi}$$, for $$h = 1, \dots, m$$, with $$\Pi_{h=1}^m x_{hi} = 1$$, which enforces the co-location requirement.

The corresponding statistic is then: $CLC_i = \Pi_{h=1}^m x_{hi} \sum_j w_{ij} \Pi_{h=1}^m x_{hj}.$

The implementation of a conditional permutation strategy follows as a direct generalization of the bivariate co-location cluster. However, for a large number of variables, such co-locations become less and less likely, and a different conceptual framework may be more appropriate.

#### Quantile Local Spatial Autocorrelation

The local Join Count statistics are designed to deal with binary variables. However, as we illustrate with our example, they can also be applied to assess the spatial autocorrelation between quantiles of a continuous variable, something that can be referred to as quantile local spatial autocorrelation. As mentioned, the continuous (linear) association between two variables that is measured by the bivariate local Moran suffers from the problem of in-situ correlation. The quantile local spatial autocorrelation sidesteps this problem by converting the continuous variable to a binary variable that takes the value of 1 for a specific quantile. For example, as in our illustration, the spatial assocation between two continuous variables can be simplified by focusing on the local autocorrelation between binary variables for the upper quintile. Since the bivariate (or multivariate) local Join Count statistic enforces co-location, the problem of in-situ correlation is controlled for. This provides an alternative to other bivariate and multivariate local spatial autocorrelation statistics, albeit with a loss of information (going from the continuous variable to a binary category).

#### Implementation

The bivariate and multivariate local Join Count functionality for the co-location case is invoked from the Cluster Maps toolbar icon, or from the menu as Space > Multivariate Local Join Count, as shown in Figure 24. Figure 24: Co-Location bivariate and multivariate local Join Count from the toolbar

The variable selection dialog takes on a special form, to allow multiple variables to be chosen. In our example, we select infq5 (top quintile for Infants) and donq5 (top quintile for Donations) as the two variables, with the spatial weights set to guerry_85_q, as in Figure 25. Figure 25: Co-Location bivariate local Join Count variable selection

Again, there is only the significance map as the result option, which gives the map in Figure 26 for the default settings of 999 permutations and p = 0.05. All the usual options for local spatial autocorrelation statistics can be set as well. Figure 26: Co-Location bivariate local Join Count significance map

The results show three significant locations. With the Save Results option, we can again add the join counts, number of neighbors and pseudo p-value to the data table. This is illustrated in Figure 27 for the three significant locations in our example. Figure 27: Co-Location bivariate local Join Count significant locations

In each instance, only two of the neighbors belong to the same category, but given the small number of neighbors, that is sufficient to indicate (weak) significance. The locations that were identified as significant constitute clusters where high values (in the upper quintile) for both variables co-locate.

## Multivariate Local Geary

In Anselin (2018), a multivariate extension of the local Geary statistic is proposed. This statistic measures the extent to which neighbors in multiattribute space (i.e., data points that are close together in the multidimensional variable space) are also neighbors in geographical space. While the mathematical formalism is easily extended to many variables, in practice one quickly runs into the curse of dimensionality.

In essence, the multivariate local Geary statistic measures the extent to which the average distance in attribute space between the values at a location and the values at its neighboring locations are smaller or larger than what they would be under spatial randomness. The former case corresponds to positive spatial autocorrelation, the latter to negative spatial autocorrelation.

An important aspect of the multivariate statistic is that it is not simply the superposition of univariate statistics. In other words, even though a location may be identified as a cluster using the univariate local Geary for each of the variables separately, this does not mean that it is also a multivariate cluster. The univariate statistics deal with distances in attribute space projected onto a single dimension, whereas the multivariate statistics are based on distances in a higher dimensional space. The multivariate statistic thus provides an additional perspective to measuring the tension between attribute similarity and locational similarity.

### Preliminaries

To illustrate the multivariate local Geary statistic, we will continue to use a bivariate example with the variables Infants (Infants) and Donations (Donatns) from the Guerry sample data set. A general sense of their overall spatial distribution can be gleaned from Figures 21 and 22 (ignoring the selection).

To provide some further perspective, the cluster maps from a univariate local Geary statistic are provided in Figures 28 and 29. The weights matrix is queen contiguity (guerry_85_q), and all the standard defaults are used, i.e., 999 permutations and p = 0.05. Figure 28: Univariate local Geary cluster map - Infants Figure 29: Univariate local Geary cluster map - Donations

### Concept

The multivariate local Geary statistic is formally the sum of individual local Geary statistics for each of the variables under consideration. For example, with $$m$$ variables, indexed by $$h$$, the corresponding expression is: $c_i = \sum_{h=1}^m \sum_j w_{ij} (x_{hi} - x_{hj})^2.$ This measure corresponds to a weighted average of the squared distances in multidimensional attribute space between the values observed at a given geographic location $$i$$ and those at its geographic neighbors $$j$$.

Inference can again be based on a conditional permutation approach. This consists of holding the tuple of values observed at $$i$$ fixed, and computing the statistic for a number of permutations of the remaining tuples over the other locations. This results in an empirical reference distribution that represents a computational approach at obtaining the distribution of the statistic under the null. The resulting pseudo p-value corresponds to the fraction of statistics in the empirical reference distribution that are equal to or more extreme than the observed statistic.

Such an approach suffers from the same problem of multiple comparisons mentioned for all the other local statistics. In addition, there is a further complication. When comparing the results for $$m$$ univariate Local Geary statistics, these multiple comparisons need to be accounted for as well. For example, for each univariate test, the target p-value of $$\alpha$$ would typically be adjusted to $$\alpha / m$$ (with $$m$$ variables, each with a univariate test), as a Bonferroni bound. Since the multivariate statistic is in essence a sum of the statistics for the univariate cases, this would suggest a similar approach by dividing the target p-value by the number of variables ($$m$$). Alternatively, and preferable, a FDR strategy can be pursued. The extent to which this actually compensates for the two dimensions of multiple comparison (multiple variables and multiple observations) remains to be further investigated.8

### Implementation

The multivariate local Geary functionality is invoked from the Cluster Maps toolbar icon, or from the menu as Space > Multivariate Local Geary, as shown in Figure 30. Figure 30: Multivariate local Geary from the toolbar

The variable selection dialog takes on the same form as for the multivariate local Join Counts, to allow multiple variables to be chosen. In our example, we select Donatns (Donations) and Infants as the two variables, with the spatial weights set to guerry_85_q, as in Figure 31. Figure 31: Multivariate local Geary variable selection

A dialog offers the option of a significance map and a cluster map, with the latter as the default. To illustrate the two types of results, we check both options.

With all the other options set to their defaults (999 permutations, p = 0.05), this yields the significance map shown in Figure 32 and the cluster map in Figure 33. Figure 32: Multivariate local Geary significance map Figure 33: Multivariate local Geary cluster map

Note that for this p-value, we obtain 24 positive clusters as well as 2 negative ones. For the latter, the distance in multi-attribute space to the neighboring locations is larger than would be the case under randomness.

All the usual options regarding the number of permutations, the significance filter, saving the results, etc. hold as for any other local spatial statistic in GeoDa. While we do not pursue this further here, it should be noted that a sensitivity analysis with respect to different significance levels is highly recommended.

### Interpretation

The results of the multivariate local Geary significance or cluster map are not always easy to interpret. Mainly, this is because we tend to simply superimpose the results of the univariate tests, whereas the multivariate statistic involves different tradeoffs.

As mentioned, the problem of multiple comparisons means that the use of the default p-value of 0.05 is almost certainly misleading. Instead, it is preferred to use the False Discovery Rate, which corresponds to about 0.0065. For example, this yields the slightly more conservative result in Figure 34, with 11 locations deemed significant. Figure 34: Multivariate local Geary significance map - False Discovery Rate

Comparing this to the two univariate local Geary cluster maps, also using the FDR criterion, reveals some interesting patterns, shown in Figure 35.

In the map for Donations, with only 7 significant locations, only three of those locations are shared with the bivariate map. Also, in the map for Infants, with only 2 significant locations, only one of those (the spatial outlier) is shared with the bivariate map. In other words, the bivariate cluster map seems to be more lenient in terms of identifying significant locations and provides additional insights beyond the univariate maps.9 Figure 35: Multivariate compared to Univariate local Geary

1. University of Chicago, Center for Spatial Data Science –

2. To recap, $$\beta = \sum_i O_i / \sum_i P_i$$, where $$O_i$$ is the number of events at $$i$$ and $$P_i$$ is the population at risk. The estimate of $$\alpha = [\sum_i P_i ( r_i - \beta )^2 ] / P - \beta / ( P / n),$$, with $$n$$ as the total number of observations, such that $$P/n$$ is the average population. Note that the estimate of $$\alpha$$ can be negative, in which case it is set to zero.

3. Load the data set into the Connect to Data Source interface in the usual way.

4. Use Map > Unique Values Map, and specify popplus as the variable. To get the same map as in the Figure, change the color of the zero category to white (right click on the legend box and set the Color for Category option to white).

5. Since the permutation test is designed to identify rare occurrences, it will mistakenly flag as significant locations that have fewer neighbors of the predominant type than expected randomly. This is because the proportion of such values in the population is larger than 0.5 so that it is likely to have a large number of neighbors of the type.

6. If you need to create this weights file, use the weights manager with NID as the ID variable.

7. Recall that the variables in the Guerry data set are scaled such that a larger value is better. In other words, large values for Infants correspond to a low degree of children born out of wedlock.

8. See Anselin (2018) for further technical discussion.

9. Note that when a significance value of p=0.05 is used for the univariate map, each of the 11 locations in the bivariate map is also significant for one or both of the univariate maps, so similar information is being captured.