# Chapter 11 Global Spatial Autocorrelation 2

## Introduction

This notebook cover the functionality of the Global Spatial Autocorrelation 2 section of the GeoDa workbook. We refer to that document for details on the methodology, references, etc. The goal of these notes is to approximate as closely as possible the operations carried out using GeoDa by means of a range of R packages.

The notes are written with R beginners in mind, more seasoned R users can probably skip most of the comments on data structures and other R particulars. Also, as always in R, there are typically several ways to achieve a specific objective, so what is shown here is just one way that works, but there often are others (that may even be more elegant, work faster, or scale better).

For this notebook, we use Cleveland house price data. Our goal in this lab is show how to assign spatial weights based on different distance functions.

### 11.0.1 Objectives

After completing the notebook, you should know how to carry out the following tasks:

Visualize bivariate spatial correlation

Assess different aspects of bivariate spatial correlation

Visualize bivariate spatial correlation among time-differenced values

Correct Moran’s I for the variance instability in rates

#### 11.0.1.1 R Packages used

**sf**: To read in the shapefile and make queen contiguity weights**spdep**: To create spatial weights structure from neighbors structure.**ggplot2**: To make customized plots such as a bivariate Moran’s I scatter plot**robustHD**: To compute standarized scores for variables and lag variables. in construction of a Moran’s I scatterplot

#### 11.0.1.2 R Commands used

Below follows a list of the commands used in this notebook. For further details and a comprehensive list of options, please consult the R documentation.

**Base R**:`install.packages`

,`library`

,`setwd`

,`summary`

,`attributes`

,`lapply`

,`class`

,`length`

,`lm`

**sf**:`st_read`

,`st_relate`

,`st_crs`

,`st_transform`

**spdep**:`nb2listw`

,`lag.listw`

**ggplot2**:`ggplot`

,`geom_smooth`

,`geom_point`

,`xlim`

,`ylim`

,`geom_hline`

,`geom_vline`

,`ggtitle`

**robustHD**:`standardized`

## 11.1 Preliminaries

Before starting, make sure to have the latest version of R and of packages that are compiled for the matching version of R (this document was created using R 3.5.1 of 2018-07-02). Also, optionally, set a working directory, even though we will not
actually be saving any files.^{29}

### 11.1.1 Load packages

First, we load all the required packages using the `library`

command. If you don’t have some of these in your system, make sure to install them first as well as
their dependencies.^{30} You will get an error message if something is missing. If needed, just install the missing piece and everything will work after that.

```
library(sf)
library(spdep)
library(ggplot2)
library(robustHD)
```

### Obtaining the Data from the GeoDa website

To get the data for this notebook, you will and to go to US County Homocide Data The download format is a
zipfile, so you will need to unzip it by double clicking on the file in your file
finder. From there move the resulting folder titled: nyc into your working directory
to continue. Once that is done, you can use the **sf** function: `st_read()`

to read
the shapefile into your R environment.

`<- st_read("natregimes/natregimes.shp") counties `

```
## Reading layer `natregimes' from data source
## `/Users/runner/work/handsonspatialdata/handsonspatialdata/natregimes/natregimes.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3085 features and 73 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.7314 ymin: 24.95597 xmax: -66.96985 ymax: 49.37173
## Geodetic CRS: WGS 84
```

### 11.1.2 Creating the weights

The **sf** neighbor functions assume that coordinates for the data are planar. Since
The current projection is that for latitude and longitude, it is best to tranform
the coordinates to a planar projection for accuracy. To check the current projection,
we use `st_crs`

. To transform the projection, we use `st_transform`

and the epsg
number of our desired projection. For more indepth information on this, please check
out the Spatial Data Handling notebook.

`st_crs(counties)`

```
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
```

`<- st_transform(counties, crs = "ESRI:102003") counties `

To start we create a function for queen contiguity, which is just `st_relate`

with
the specified pattern for queen contiguity which is `F***T****`

`<- function(a, b = a) st_relate(a, b, pattern = "F***T****") st_queen `

We apply the queen contiguity function to the voronoi polygons and see that the class
of the output is **sgbp**. This structure is close to the **nb** structure, but has
a few difference that we will need to correct to use the rest of **spdep** functionality.

```
<- st_queen(counties)
queen.sgbp class(queen.sgbp)
```

`## [1] "sgbp" "list"`

This function converts type **sgbp** to **nb**. It is covered in more depth in the
Contiguity Based Weight notebook. In short, it explicitly changes the name of the
class and deals with the observations that have no neighbors.

```
<- function(x, ...) {
as.nb.sgbp <- attributes(x)
attrs <- lapply(x, function(i) { if(length(i) == 0L) 0L else i } )
x attributes(x) <- attrs
class(x) <- "nb"
x }
```

`<- as.nb.sgbp(queen.sgbp) queen.nb `

To go from neighbors object to weights object, we use `nb2listw`

, with default parameters, we will
get row standardized weights.

`<- nb2listw(queen.nb) queen.weights `

## 11.2 Bivariate spatial correlation - a Word of Caution

The concept of bivariate spatial correlation is complex and often misinterpreted. It is typically considered to the correlation between one variable and the spatial lag of another variable, as originally implemented in the precursor of GeoDa. However this does not take into account the inherent correlation between the two variables. More precisely, the bivariate spatial correlation is between \(x_i\) and \(\Sigma_jw_{ij}y_j\), but does not take into account the correlation between \(x_i\) and \(y_i\), i.e. between two variables at the same location.

As a result, this statistic is often interpreted incorrectly, as it may overestimate the spatial aspect of the correlation that instead may be due mostly to in-place correlation.

Below, we provide a more in-depth assessment of the different aspects of bivariate spatial and non-spatial association, but first we turn to the original concept of a bivariate Moran scatter plot.

## 11.3 Bivariate Moran Scatter Plot

### 11.3.1 Concept

In its initial conceptualization, as mentioned above, a bivariate Moran scatter plot extends the idea of a Moran scatter plot with a variable on the x-axis and its spatial lag on a y-axis to a bivariate context. The fundamental difference is that in the bivariate case the spatial lag pertains to a different variable. In essence, this notion of bivariate spatial correlation measures the degree to which the value for a given variable at a location is correlated with its neighbors for a different variable.

As in the univariate Moran scatter plot, the interest is in the slope of the linear fit. This yields a Moran’s I-like statistic as:

\[I_B=\frac{\Sigma_i(\Sigma_jw_{ij}*x_i)}{\Sigma_ix_i^2}\]

or, the slope of a regression of Wy on x.

As before, all variables are expressed in standardized form, such that their means are zero and their variance one. In addition, the spatial weights are row-standardized.

Note that, unlike in the univariate autocorrelation case, the regression of x on Wy also yields an unbiased estimate of the slope, providing an alternative perspective on bivariate spatial correlation.

A special case of bivariate spatial autocorrelation is when the variable is measured at two points in time, say \(z_{i,t}\) and \(z_{i,t-1}\). The statistic then pertains to the extent to which the value obsereved at a location at a given time, is correlated with its value at neighboring locations at a different point in time.

The natural interpretation of this concept is to relate \(z_{i,t}\) to \(\Sigma_jw_{ij}z_{j,t-1}\), i.e. the correlation between a value at t and its neighbors in a previous time period:

\[I_T=\frac{\Sigma_i(\Sigma_jw_{ij}z_{j,t-1}*z_{i,t})}{\Sigma_iz_{i,t}^2}\] i.e., the effect neighbors in t-1 on the present value.

Alternatively, and maybe less intuitively, one can relate the value at a previous time period \(z_{t-1}\) to its neighbors in the future, \(\Sigma_jw_{ij}z_t\) as:

\[I_T=\frac{\Sigma_i(\Sigma_jw_{ij}z_{j,t}*z_{i,t-1)}}{\Sigma_iz_{i,t-1}^2}\] i.e., the effect of a location at t-1 on its neighbors in the future.

While formally correct, this may not be a proper interpretation of the dynamics involved. In fact, the notion of spatial correlation pertains to the effect of neighbors on a central location, not the other way around. While the Moran scatter plot seems to reverse this logic, this is purely a formalism, without any consequences in the univariate case. However, when relating the slope in the scatter plot to the dynamics of a process, this should be interpreted with caution.

A possible source of confusion is that the proper regression specification for a dynamic process would be as:

\[z_{i,t}=\beta_1\Sigma_jw_{ij}z_{j,t-1} +u_i\] with \(u_i\) as the usual error term, and not as:

\[z_{i,t-1}= \beta_2\Sigma_jw_{ij}z_{j,t} + u_i\]

which would have the future predicting the past. This contrasts with the linear regression specification used (purely formally) to estimate the bivariate Moran’s I, for example:

\[\Sigma_jw_{ij}z_{j,t-1} = \beta_3z_{i,t} + u_i\]

In terms of the interpretation of a dynamic process, only \(\beta_1\) has intuitive appeal. However, in terms of measuring the degree of spatial correlation between past neighbors and a current value, as measured by a Moran’s I coefficient, \(\beta_3\) is the correct interpretation.

As mentioned earlier, a fundamental difference with the univariate case is that the spatially lagged variable is no longer endogenous, so both specifications can be estimated by means of ordinary least squares.

In the univariate case, only the specification with the spatially lagged variable on the left hand side yields a valid estimate. As a result, for the univariate Moran’s I, there is no ambiguity about which variables should be on the x-axis and y-axis.

Inference is again based on a permutation approach, but with an important difference. Since the interest focuses on the bivariate spatial association, the values for x are fixed at their locations, and only the values for y are randomly permuted. In the usual manner, this yields a reference distribution for the statistic under the null hypothesis that the spatial arrangement of the y values is random. As mentioned, it is important to keep in mind that since the focus is on the correlation between the x value at i and the y values at the neighboring locations, the correlation between x and y at location i is ignored.

### 11.3.2 Creating a bivariate Moran scatter plot

To create the bivariate Moran scatterplot, we will need to standardize the variables beofre proceeding. This
is simply, putting the variables in terms of z scores. To do this, we use `standardize`

from the
**robustHD** package.

```
$hr80st <- standardize(counties$HR80)
counties$hr90st <- standardize(counties$HR90) counties
```

Now that we have standardized the variables, we can create the lag variable for homocide rates in 1980.
This is done with `lag.listw`

from **spdep**. For more a more indepth coverage of this function , please see
the Application of Spatial Weights notebook.

`$whr80st <- lag.listw(queen.weights,counties$hr80st) counties`

In constructing the majority of our plots, we use **ggplot2** because it makes aesthically pleasing and functional
plots. The purpose of this tutorial is not to teach **ggplot2**, but to use it in creqating specializzed plots
such as the bivariate Moran plot. For indepth coverage of basic **ggplot2** functionality please see the
Exploratory Data Analysis notebooks or ggplot2 documentation.

To make the bivariate Moran, we make a scattplot of the lag variable for 1980 homocide rates vs the homocide
rate for 1990. In GeoDa the lag variable and standardization are done automatically and all you have to do is
select the desired variables. The custom specifications use for this plot are the dashed lines at the x and y
axis and the simple regression line. `geom_vline`

and `geom_hline`

.

```
ggplot(data = counties, aes(x=hr90st,y=whr80st)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty =2) +
xlim(-12,12) +
ylim(-12,12) +
ggtitle("Bivariate Moran's I, hr(90) and lagged hr(80) ")
```

To get the summary statistics and bivariate Moran’s I statistic, we will run a rgeression separate from the plot.
For this regression, we use the same variables has the bivariate Moran plot in the base R `lm`

function.
To view the results, we use the `summary`

command. The Moran’s I statistic is the slope of the regression
line, which is the coefficient of **hr90st** in the `summary`

results

```
<- lm(whr80st ~ hr90st, data = counties)
lmfit summary(lmfit)
```

```
##
## Call:
## lm(formula = whr80st ~ hr90st, data = counties)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9290 -0.4171 -0.1024 0.3595 2.6279
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.004734 0.010694 0.443 0.658
## hr90st 0.359668 0.010696 33.628 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.594 on 3083 degrees of freedom
## Multiple R-squared: 0.2684, Adjusted R-squared: 0.2681
## F-statistic: 1131 on 1 and 3083 DF, p-value: < 2.2e-16
```

### 11.3.3 A closer look at bivariate spatial correlation

In contrast to the univariate Moran scatter plot, where the interpretation of the linear fit is unequivocably Moran’s I, there is no such clarity in the bivariate case. In addition to the interpretation offered above, which is a traditional Moran’s I-like coefficient, there are at least four additional perspectives that are relevant. We consider each in turn.

#### 11.3.3.1 Serial (temporal) correlation

To examine the temporal correlation of homocide rates, we plot the standardized variable for each time period. This gives us a sense of the temporal autocorrelation of homocide rates for a lag order of 1.

```
ggplot(data = counties, aes(x = hr80st, y = hr90st)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty =2)
```

Again we us the `lm`

function and `summary`

function to get a sense of the relationship
between the variables. With a highly significant slope of .552, strong temporal correlation
is suggested.

```
<- lm(hr90st ~ hr80st, data = counties)
lmfit summary(lmfit)
```

```
##
## Call:
## lm(formula = hr90st ~ hr80st, data = counties)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8479 -0.3701 -0.1609 0.2806 10.3772
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.018e-16 1.501e-02 0.00 1
## hr80st 5.526e-01 1.501e-02 36.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8336 on 3083 degrees of freedom
## Multiple R-squared: 0.3054, Adjusted R-squared: 0.3052
## F-statistic: 1356 on 1 and 3083 DF, p-value: < 2.2e-16
```

#### 11.3.3.2 Serial (temporal) correlation between spatial lags

Next we run a regression of the spatial lags of homocide rates for both time periods. This will be a regression of \(\Sigma_jw_{ij}y_{j,t}\) on \(\Sigma_jw_{ij}y_{j,t-1}\). To begin, wewill need to create the lag variable for 1990.

`$whr90st <- lag.listw(queen.weights, counties$hr90st) counties`

To make this plot, we just plot the lag variable for 1990 against the one for 1980.

```
ggplot(data = counties, aes(x= whr80st, y = whr90st)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty =2)
```

Again the slope is highly signiificcant with a positive slope of .801, even strong than the actual variable’s temporal correlation.

```
<- lm(whr90st ~ whr80st, data = counties)
lmfit summary(lmfit)
```

```
##
## Call:
## lm(formula = whr90st ~ whr80st, data = counties)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.73241 -0.20479 -0.03393 0.16637 2.13011
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0005499 0.0074544 -0.074 0.941
## whr80st 0.8014254 0.0107383 74.633 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.414 on 3083 degrees of freedom
## Multiple R-squared: 0.6437, Adjusted R-squared: 0.6436
## F-statistic: 5570 on 1 and 3083 DF, p-value: < 2.2e-16
```

#### 11.3.3.3 Space-time regression

In this case we have already computed the variables need to construct the plot, so we
can start with **ggplot2**. The regression here is \(y_{i,t}\) on \(\Sigma_jw_{ij}y_{j,t-1}\)

```
ggplot(data = counties, aes(x= whr80st, y = hr90st)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty =2)
```

Now we run the `lm`

function to get the regression statistics.

```
<- lm(hr90st ~ whr80st, data = counties)
lmfit summary(lmfit)
```

```
##
## Call:
## lm(formula = hr90st ~ whr80st, data = counties)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3519 -0.4569 -0.1702 0.2711 10.1136
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.003532 0.015403 -0.229 0.819
## whr80st 0.746141 0.022188 33.628 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8555 on 3083 degrees of freedom
## Multiple R-squared: 0.2684, Adjusted R-squared: 0.2681
## F-statistic: 1131 on 1 and 3083 DF, p-value: < 2.2e-16
```

This yields a highly signifcant slope of .746

#### 11.3.3.4 Serial and space-time regression

Finally, we can include both in-situ correlation and space-time correlation in a
regression specification of \(y_{i,t}\) i.s **hr90st** on both \(x_{i,t-1}\) (ie, **hr80st**)
and \(\Sigma_jw_{ij}x_{j,t-1}\)(ie **whr80st**)

```
<- lm(hr90st ~ whr80st + hr80st, data = counties)
lmfit summary(lmfit)
```

```
##
## Call:
## lm(formula = hr90st ~ whr80st + hr80st, data = counties)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2802 -0.3807 -0.1202 0.2574 10.3821
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002093 0.014282 -0.147 0.883
## whr80st 0.442224 0.024628 17.956 <2e-16 ***
## hr80st 0.383861 0.017099 22.450 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7932 on 3082 degrees of freedom
## Multiple R-squared: 0.3712, Adjusted R-squared: 0.3708
## F-statistic: 909.7 on 2 and 3082 DF, p-value: < 2.2e-16
```

In the regression, both effects are highly significant, with a coefficient of 0.384 for the in-situ effect, and 0.442 for the space-time effect.

Needless to say, focusing the analysis solely on the bivariate Moran scatter plot provides only a limited perspective on the complexity of space-time (or, in general, bivariate) spatial and space-time associations.

## 11.4 Differential Moran Scatter Plot

### 11.4.1 Concept

An alternative to the bivariate spatial correlation between a variable at one point in time and its neighbors at a previous (or, in general, any different) point in time is to control for the locational fixed effects by diffencing the variable. More precisely, we apply the computation of the spatial autocorrelation statistic to the variable \(y_{i,t} - y_{i,t-1}\).

This takes care of any temporal correlation that may be the result of a fixed effect(i.e. a set of variables determining the value of \(y_i\) that remain unchanging over time). Specifically if \(\mu_i\) were a fixed effect associated with location i, we could express the value at each location for time t as the sum of some intinsic value and the fixed effect, $y_{i,t} = y_{i,t}^* + _i $, where \(\mu_i\) is unchanged over time (hence, fixed). The presence of \(\mu\) induces a temporal correlation between \(y_{i,t}\) and \(y_{i,t-1}\), above and beyond the correlation between \(y_{i,t}^*\) and \(y_{i,t-1}^*\). Taking the difference eliminates the fixed effect and ensures that any remainin correlation is solely due to \(y^*\).

A differential Moran’s I is then the slope in a regression of the spatial lag of the difference i.e. \(\Sigma_jw_{ij}(y_{j,t} - y_{j,t-1})\) on the difference \(y_{i,t} - y_{i,t-1}\).

All the caveats and different interpretations that apply to the generic bivariate case are also relevant here.

### 11.4.2 Implementation

Same as the other sections, we must first calculate and standardize the variables
before proceeding to construct the plot. In GeoDa this process is done automatically,
but here we will need a few extra steps. First we calcualte the standardized score of
the difference between the variable, NOT the difference between the standardized
score of each variable. To do this, we just use `standardize`

on the difference
between **HR80** and **HR90**. Then we calcualte the lag variable with `lag.listw`

on the standardized difference variable.

```
$hr_diff_st <- standardize(counties$HR90 - counties$HR80)
counties$hr_diff_lag <- lag.listw(queen.weights, counties$hr_diff_st) counties
```

To construct the plot, we use the same functions as earlier, but just replace the variables. In this case, we use the ones calcualted above for the standardized difference then the lag of the standardized difference.

```
ggplot(data = counties, aes(x = hr_diff_st, y = hr_diff_lag)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty =2)
```

```
<- lm(hr_diff_lag~hr_diff_st, data = counties)
lmfit summary(lmfit)
```

```
##
## Call:
## lm(formula = hr_diff_lag ~ hr_diff_st, data = counties)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.23658 -0.21238 0.03382 0.22281 2.45484
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001690 0.008258 -0.205 0.838
## hr_diff_st 0.052007 0.008260 6.297 3.47e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4587 on 3083 degrees of freedom
## Multiple R-squared: 0.0127, Adjusted R-squared: 0.01238
## F-statistic: 39.65 on 1 and 3083 DF, p-value: 3.474e-10
```

The slope is about .052, which is much smaller than in magnitude than the other measures used thus far.

## 11.5 Moran Scatter Plot for EB Rates

### 11.5.1 Concept

An Empirical Bayes (EB) standardization was suggested by Assunção and Reis (1999) as a means to correct Moran’s I spatial autocorrelation test statistic for varying population densities across observational units, when the variable of interest is a rate or proportion. This standardization borrows ideas from the Bayesian shrinkage estimator outlined in discussion of Empirical Bayes smoothing.

This approach is different from EB smoothing in that the spatial autocorrelation is not computed for a smoothed version of the original rate, but for a transformed standardized random variable. In other words, the crude rate is turned into a new variable that has a mean of zero and unit variance, thus avoiding problems with variance instability. The mean and variance used in the transformation are computed for each individual observation, thereby properly accounting for the instability in variance.

The technical aspects are given in detail in Assunção and Reis (1999) and in the review by Anselin, Lozano-Gracia, and Koschinky (2006), but we briefly touch upon some salient issues here. The point of departure is the crude rate, say \(r_i = O_i/P_i\), where \(O_i\) is the count of events at location i and \(P_i\) is the corresponding population at risk.

The rationale behind the Assuncao-Reis approach is to standardize each \(r_i\) as

\[z_i = \frac{r_i - \beta}{\sqrt{\alpha + (\beta/P_i)}}\]

using an estimated mean \(\beta\) and standard error \(\sqrt{\alpha + (\beta/P_i)}\). The parameter \(\alpha\) and \(\beta\) are related to the prior distribution for the risk, which is detailed in the notebook on rate mapping.

In practice, the prior parameters are estimated by means of the so-called method of moments (e.g., Marshall 1991), yielding the following expressions:

\[\beta=O/P\]

With O as the total event count \((\Sigma_iO_i)\), and P as the total population count \((\Sigma_iP_i)\), and

\[\alpha = [\Sigma_iP_i(r_i - \beta)^2]/P - \beta/(P/n)\]

with n as the total number of observations (in other words, P/n is the average population)

One problem with the method of moments estimator is that the expression for \(\alpha\) could yield a negative result. In that case, its value is typically set to zero, i.e. \(\alpha = 0\). However, Assuncao and Reis(1999), the value for \(\alpha\) is set to zero when the resulting estimate for variance is negative, that is, \(\alpha + \beta/P_i <0\). Slight differences in the standardized variates may result depending on the convention used. In GeoDa, when the variance estimate is negative, the original raw rate is used.

### 11.5.2 Implementation

To calculate the empirical bayes rates, we will use basic R operations. We start by calculating the initial variables.

```
<- sum(counties$HC60) / sum(counties$PO60)
beta <- counties$HC60 / counties$PO60
r_i <- counties$PO60
P_i <- length(r_i)
n <- sum(counties$PO60) P
```

With the parameters from above, we can calculate \(\alpha\), by the formula specified in the concept section.

`<- sum(P_i * (r_i - beta)^2) / P - beta/(P/n) alpha `

With \(\alpha\), we can get the standard error used to later calculate the standardized values.

`<- sqrt(alpha + beta/P_i) se `

Now the standardized values are just the crude rate minus \(\beta\) and divided by the standard error.

`<- (r_i - beta) / se z_i `

With our standardized values, we can make the lag variable to constuct our
empirical bayes Moran plot. We use `lag.listw`

as before to make the spatial lag
varible.

```
$hc60_po60 <- z_i
counties$lag_hc60_po60 <- lag.listw(queen.weights, z_i) counties
```

Now with the same **ggplot2** functions as earlier we plot the smoothed rate
variable against the lag smoothed rate variable.

```
ggplot(data = counties, aes(x = hc60_po60, y = lag_hc60_po60)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty =2)
```

To get the regression statistics, we run `lm`

and then examine the summary with
`summary`

.

```
<- lm(lag_hc60_po60 ~ hc60_po60, data = counties)
lmfit summary(lmfit)
```

```
##
## Call:
## lm(formula = lag_hc60_po60 ~ hc60_po60, data = counties)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.05433 -0.34093 -0.09708 0.27749 2.89081
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001913 0.009724 -0.197 0.844
## hc60_po60 0.518303 0.010598 48.906 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.54 on 3083 degrees of freedom
## Multiple R-squared: 0.4369, Adjusted R-squared: 0.4367
## F-statistic: 2392 on 1 and 3083 DF, p-value: < 2.2e-16
```