Guangdong-Hong Kong-Macau Public Administration Data Analysis Contest Workshop

An Introduction to Geospatial Data Analysis

BAO Le

City University of Hong Kong

November 4, 2024

Introduction & Plan

宝乐 BAO, Le
- Assistant Professor of Political Science
- Department of Public and International Affairs
- Political Methodology, Political Behavior, & Environmental Politics

Today:
- The basics of spatial data
- Spatial autocorrelation
- Spatial interpolation
- Spatial regression

Spatial Analysis

Space and geography play prominent roles in many social science theories.
Many social science theories argue that nearby units (people, neighborhoods, or countries) are likely to behave similarly.
Galton’s Problem (1888): whether similar behaviors across regions arise due to diffusion (interaction) or because similar characteristics are located together geographically (clustering).
Why spatial?
- Traditional models use dummy variables to represent geographical differences, but this fails to capture the complex nature of spatial relationships.

空间和地理位置在很多社会科学理论中都扮演着重要角色。在社会科学领域，地理上的接近性常常被认为对人类行为有显著的影响。例如，在人口学、社会学、政治学等领域，理论通常假设相邻的单位（如个人、社区、甚至国家）会表现出相似的行为或特征。很多社会现象不是随机分布的，而是根据地理位置形成某种模式。邻近的区域通常会在行为或特征上表现出一致性，比如相邻的社区可能在经济状况、犯罪率、健康状况等方面表现出相似的特征。空间分析帮助我们揭示这些模式及其背后的原因，理解地理接近如何影响社会行为。

1888年，科学家弗朗西斯·加尔顿（Sir Francis Galton）提出了一个影响深远的问题，后来被称为“加尔顿问题”。加尔顿质疑：不同地区表现出相似行为，究竟是因为这些行为通过互动在邻近地区之间扩散，还是因为相似的特征在地理上聚集在一起？例如，在相邻的社区中，支持某个政党的比例可能相似。那么，这种相似性是因为居民通过讨论和互动影响了彼此的政治倾向，还是因为这些社区本身就聚集了相似收入水平、教育程度的人群？加尔顿问题提醒我们，在空间分析中必须区分“扩散”与“聚集”两种不同的空间依赖来源，才能更准确地理解和解释数据中的空间依赖现象。

传统的统计模型通常使用哑变量（dummy variables）来代表地理差异，比如通过不同的地区代码来区分区域间的差异。但是，这种做法有其局限性，因为它无法反映空间关系的复杂性。空间关系不是简单的地理单位的差异，而是连续性和相互影响的表现。哑变量无法捕捉邻近地区之间相互依赖的程度，也难以准确反映空间扩散或聚集的过程。空间分析为我们提供了更复杂的工具，使我们能够明确建模空间上的相关性，分析社会现象中相邻单位间的互动和依赖。这不仅增强了理论的解释力，还帮助我们识别地理因素在社会现象中的影响力，为政策制定提供了更有针对性的依据。

Spatial Dependence

Tobler’s first law of geography

“Everything is related to everything else, but near things are more related than distant things.”

                                                                              - Tobler (1970)

Examples of spatial dependence: crime patterns, disease spread, housing market, political orientations, transportation, environmental equity, etc.
We oftentimes even use spatial dependence without realizing it:

在空间分析中，有一个非常重要的概念叫做空间依赖性（spatial dependence）。简单来说，空间依赖性指的是地理上相邻的区域之间往往表现出更强的相似性。这种依赖性对于理解社会科学中的许多现象至关重要。在这里，我们可以引用托布勒的地理学第一定律来更直观地理解这个概念。

美国著名的地理学家，特别以其在地理信息科学（GIS）和空间分析方面的贡献而闻名。他在1970年提出了地理学第一定律，这一理论后来成为空间分析的基础。托布勒（Tobler）在1970年提出了一句经典的描述——“一切事物都相互关联，但近的事物比远的事物关联更强。”这句话总结了空间依赖性的基本思想。换句话说，地理上接近的地方往往会在行为、特征或结果上表现出更大的相似性。这种空间相似性可以在许多社会现象中观察到。

空间依赖性不仅是一个理论概念，它在实际应用中非常普遍。比如：

犯罪模式：在城市中，犯罪通常在特定区域内聚集形成高犯罪率的“热点”。疾病传播：传染病在地理上相邻的社区之间迅速传播，例如COVID-19的社区传播。房地产市场：房价在相邻区域之间通常呈现出相似的趋势，因为相邻的社区受到相似的经济和社会因素影响。政治倾向：选民的政治倾向往往在邻近区域中表现出一致性，因为相似的文化和社交网络会强化同一地区居民的政治观念。交通系统：区域之间的交通模式显示出空间依赖性，比如城市中的通勤路径往往呈现出空间规律。环境公平性：一些社区往往更容易受到环境污染的影响，而这些区域通常聚集了特定的社会经济群体，显示出空间上的不平等。

实际上，即使我们没有明确意识到，很多时候我们也会使用空间依赖性的信息。例如，右边这张图片是来自一项发表于《美国国家科学院院刊》（PNAS）的研究。研究人员通过深度学习和Google街景数据，分析社区的地理分布和特征，以估算出美国各地的社区人口组成。这项研究利用了空间依赖性，通过邻近性的信息来推测出每个区域的特征，从而帮助我们了解社区的社会和经济状况。

空间依赖性帮助我们在空间维度上分析和理解社会现象。无论是直接研究邻近区域的相似性，还是无意识地应用邻近性概念，空间依赖性都是社会科学中分析地理数据的关键因素。这一概念不仅让我们对数据有更深入的理解，也为我们提供了新的分析视角，使我们能够更精确地进行预测和制定政策。

Spatial Dependence

Diffusion-Based Dependence: Behaviors or phenomena in one area affect neighboring areas through direct interaction.
- E.g. Political discussions across adjoining neighborhoods can increase support for a candidate in nearby areas.
Attributional Dependence: Geographical clustering occurs because similar characteristics are grouped together.
- E.g. Adjoining areas may favor the same candidate not due to interaction, but because of shared demographic factors like income level.
Solutions:
- Spatial Lag Models: Capture direct influences from neighboring units, modeling how a change in one area directly impacts nearby areas.
- Spatial Error Models: Address spatial correlation in the error terms, often due to unobserved spatial factors that create similar outcomes across regions.

基于扩散的依赖性 (Diffusion-Based Dependence) 在环境污染治理政策中，基于扩散的依赖性非常明显，因为污染通常不会局限在单一区域，而是会扩散到周边地区。这种情况下，一个区域的污染治理效果可能会通过地理位置的邻近性影响到相邻的地区。

例子：假设一个城市的某个区域采取了严格的空气污染治理措施，如限制工业排放、推广清洁能源等。这一政策的实施可能会改善该地区的空气质量，同时也会对相邻地区的空气质量产生积极影响。邻近区域的空气污染减少并不是因为这些区域本身采取了措施，而是因为治理效果的扩散作用。这种现象便是基于扩散的空间依赖性。在这种情况下，环境治理政策的效果不仅局限在政策实施区域，还会通过扩散效应影响周边区域。因此，政策制定者在设计环境治理政策时，应该考虑如何在更大范围内实现环保效益，从而更好地发挥政策的区域协同作用。

属性聚集的依赖性 (Attributional Dependence) 另一方面，一些区域的污染状况和治理效果可能是由于其相似的地理属性或经济结构导致的，而不是区域间的直接相互影响。此类依赖性称为属性聚集的依赖性。

例子：在一些工业区集中的地区，空气污染水平往往较高。这些区域之间可能并没有直接的互动影响，而是因为它们都拥有相似的产业结构（如制造业密集），导致空气污染的程度相似。即使没有污染治理政策，这些相似的特征也会使得这些区域的环境状况相近。这种属性聚集的依赖性提醒政策制定者，不同区域的环境治理需要因地制宜。例如，在多个高污染工业区集中分布的城市，政策制定者可能需要制定特定的环保法规和监管措施来解决这些地区的共性问题，而不仅仅依赖一刀切的政策措施。

解决方法：空间滞后模型与空间误差模型为了更准确地分析环境政策的空间依赖性，我们可以使用空间滞后模型和空间误差模型。

空间滞后模型 (Spatial Lag Models)：空间滞后模型能够帮助我们量化一个区域的治理政策对周边区域的直接影响。例如，通过模型可以分析一个区域的空气污染控制政策如何在空间上扩散，并影响周边区域的空气质量。

空间误差模型 (Spatial Error Models)：当污染状况的相似性主要由区域特征的相似性引起时，空间误差模型可以帮助我们控制这些因素带来的误差。这样，模型就能更准确地评估治理政策的效果，而不会受到未观察到的因素（如产业结构相似性）的干扰。

总结在环境污染治理中，理解空间依赖性能够帮助我们更有效地设计政策。无论是考虑扩散效应还是区域间的共性特征，空间统计模型都为我们提供了有力的分析工具，使得政策制定能够更加精准，有助于实现更广泛的环境改善目标。

Basics: Types of Spatial Data

Vector data
- Spatial objects are often depicted using vector data, which details the object’s “geometry” or “shape” and includes additional variables, or “attributes”.
- At its simplest level, vector data comprises of individual points stored as coordinate pairs that indicate a physical location in the world.

name	longitude	latitude	aqi
River Terrace	-76.95807	38.89557	42
Rockville	-77.10688	39.11431	30
Springfield	-77.18347	38.76835	34

在空间分析中，我们经常会处理不同类型的空间数据。今天，我们首先介绍矢量数据（Vector Data），这是一种在地理信息系统（GIS）中非常常见的数据类型。矢量数据帮助我们表示空间对象的几何形状和位置，是我们进行地理分析的基础工具之一。

什么是矢量数据？矢量数据通过一系列的坐标点来表示空间对象的位置和形状。它通常包含两部分信息：

几何形状（Geometry）：即物体在地理空间中的形状或位置。简单来说，几何形状可以通过一组坐标点来表示，这些点指示出空间对象在世界中的具体位置。属性数据（Attributes）：除了几何信息，矢量数据还包含关于对象的其他特征信息，比如名称、分类、数值等。这些属性信息帮助我们更全面地理解空间对象的特性。例如，在矢量数据中，一条河流可以用一组坐标点来描绘其流向，同时还可以附加河流的长度、流量等信息作为属性。

矢量数据的基本类型在最简单的层次上，矢量数据可以由点数据（Point Data）组成。每个点由坐标对（如经度和纬度）表示，代表着空间中的一个具体位置。通过这些点，我们可以构建更复杂的空间对象，比如线（Line）和面（Polygon）：

点：表示特定位置的空间数据，比如一个城市或一个景点的坐标。线：一系列相连的点可以表示一条路径或河流等线性对象。面：由多条线围成的区域可以表示一个特定的地理范围，如国家、湖泊或城市区域。 3. 示例代码解释在右侧的代码示例中，我们用R语言创建了一个简单的数据框，存储了三个位置的坐标和空气质量指数（AQI）：

在这个数据框中，我们定义了三个地点的名称（name）、经度（longitude）、纬度（latitude）和空气质量指数（AQI）。这些数据可以作为矢量数据来进行地理分析，比如将它们绘制在地图上以查看空气质量的空间分布情况。

Vector Data

plot(dmv_air$longitude, dmv_air$latitude, pch=20, 
     col=heat.colors(3)[rank(dmv_air$aqi)], 
     xlim=c(-77.25, -76.75), main='Air Quality Index', 
     ylab = "Latitude", xlab = "Longitude", bg="white")
text(dmv_air$longitude, dmv_air$latitude, labels=dmv_air$name, 
     pos=4, cex=0.7, offset = 0.5)
legend("topright", legend=unique(aqi), fill=heat.colors(3)[order(unique(aqi))], 
     title="AQI", bg="white")

Basics: Types of Spatial Data

Raster data
- Raster data is typically used to depict continuous spatial phenomena (e.g. altitude or temperature gradients).
- A raster dataset subdivides the spatial area into a grid of uniform rectangles, each holding one or more values for the variables of interest.
- In raster data, geometry is inherently encoded rather than explicitly stored as coordinates.

除了矢量数据之外，栅格数据（Raster Data）是空间数据的另一种重要类型。栅格数据通常用于表示连续的空间现象，比如海拔、温度等变量。接下来，我们将深入探讨栅格数据的特点和应用场景。

什么是栅格数据？栅格数据将空间区域划分为一系列大小均匀的矩形网格，每个网格单元（或称为像素）包含一个或多个数值，用于表示该位置的属性或状态。

连续现象：栅格数据非常适合表示连续分布的空间现象。这些现象通常没有明确的边界，而是逐渐变化的，比如温度变化、降水量分布、海拔等。通过栅格数据，我们可以创建出地理现象的“地图”或“地形”，反映变量的空间变化趋势。单一值：每个网格单元都包含一个数值，用于代表该区域的特定属性。比如，在温度栅格数据中，每个单元格可能存储当前温度值，而在海拔数据中，每个单元格则存储该位置的海拔高度。 2. 栅格数据的结构特点在栅格数据中，几何形状并不是通过坐标来明确存储的，而是隐含在网格结构中。具体来说：

每个网格单元代表空间中的一个小区域，整个栅格数据就像一个“像素化”的地图。栅格大小（或分辨率）决定了地图的精度——分辨率越高，地图越精细，但文件也会越大。栅格数据的每一行和每一列都对应地理空间中的特定位置，栅格的行列位置与地理坐标之间有固定的换算关系。这意味着我们不需要为每个单元格单独存储坐标，数据的几何信息隐含在整体结构中。 3. 栅格数据的应用场景栅格数据广泛用于地理信息系统（GIS）和环境科学中，适合处理以下类型的空间数据：

海拔：用于表示地形高度，例如山脉、平原和河谷的分布。温度分布：反映不同地理区域的温度变化，可以用于气象预测和气候分析。土壤湿度或植被指数：这些数据通常在农业、生态学和环境管理中应用，帮助监测土地健康状况和生态系统的变化。通过栅格数据，我们能够快速处理和分析大规模的空间数据，捕捉到地理现象的微小变化。例如，卫星图像通常以栅格数据形式存储，每个像素代表地球表面的一小块区域，可以用于遥感监测和地表分析。

总结栅格数据提供了一种有效的方式来表示和分析连续的空间现象。它通过将空间区域划分为均匀的网格，将每个位置的属性数据隐含在网格结构中，为我们提供了丰富的地理信息。无论是在气候研究还是生态监测中，栅格数据都发挥着至关重要的作用。

Raster Data

library(raster)
dc_extent <- extent(-77.12, -76.91, 38.79, 39.05)
ras <- raster(dc_extent, nrow=25, ncol=25)
values(ras) <- runif(ncell(ras))
ras

class      : RasterLayer 
dimensions : 25, 25, 625  (nrow, ncol, ncell)
resolution : 0.0084, 0.0104  (x, y)
extent     : -77.12, -76.91, 38.79, 39.05  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 0.0001935221, 0.9998556  (min, max)

Raster Data

plot(dc_extent, col=NA)
plot(ras, add=T, legend=F)

Vector Data

Vector data consists of points, lines, and polygons, etc.
- For example, a point might represent a location of air quality sensor with attributes including the observation date, the observed air quality, and the pollutant’s information.
- A line, a series of connected coordinate pairs, could represent a river system.
- Polygons, which are closed polylines, are often used to represent areas such as geographical areas.

在空间分析中，矢量数据是表示空间对象的一种主要形式。矢量数据的特点是由点、线和多边形等几何形状组成，能够精准地表示地理对象的位置和边界。接下来，我们将深入了解矢量数据的三种基本类型：点、线和多边形，并通过一些实际例子来说明它们的应用。

点 (Point) 在矢量数据中，点是一种最简单的空间对象。每个点由一对坐标（通常是经纬度）来表示具体的位置。点不仅可以用于表示地点，还可以附带一些属性信息，用于描述该位置的特定特征。

例子：一个空气质量传感器的位置可以用一个点来表示，同时还可以附加属性信息，例如观测日期、空气质量指数（AQI）以及特定污染物的信息。这些属性数据帮助我们更全面地了解该位置的环境状况。 2. 线 (Line) 线是由一系列连续的坐标点连接而成，用于表示具有长度的地理对象。线可以用于描述各种路径、流动或延展性的空间对象。

例子：一条河流可以用线来表示，河流的流向和路径通过连接的坐标点清晰地展示出来。除了河流，交通路线、道路和地铁线路等也常用线来表示。这种空间表示方式有助于我们分析这些地理对象的延伸方向和结构。 3. 多边形 (Polygon) 多边形是一种由闭合的多条线段构成的几何形状，用于表示有面积的空间对象。多边形在矢量数据中应用广泛，尤其适合表示具有边界的区域。

例子：例如，我们可以用多边形来表示城市、国家或某个生态保护区的地理范围。每个多边形不仅定义了该区域的边界，还可以附带区域的其他属性信息，比如人口密度、土地利用类型等。多边形数据帮助我们更直观地展示和分析区域的空间特征和属性。

Line

library(sp)
line_coords <- data.frame(
  longitude = c(-77.3, -77.2086, -77.2367, -77.2555, -77.2865, -76.8),
  latitude = c(38.92, 38.9041, 38.8954, 38.8816, 38.8688, 38.85)
)

line <- SpatialLines(list(Lines(Line(line_coords), ID="segment")))
line

class       : SpatialLines 
features    : 1 
extent      : -77.3, -76.8, 38.85, 38.92  (xmin, xmax, ymin, ymax)
crs         : NA

Line

plot(line, col="blue", lwd=2, main="Potomac River", 
     ylab="Latitude", xlab="Longitude")

Polygon

polygon_coords <- rbind(
    c(-77.11976, 38.93434),
    c(-77.04102, 38.99555),
    c(-77.00255, 38.96553),
    c(-76.90939, 38.89285),
    c(-76.97950, 38.83781),
    c(-77.03901, 38.79164),
    c(-77.03907, 38.84127),
    c(-77.03910, 38.86811),
    c(-77.09020, 38.90421),
    c(-77.11976, 38.93434)
)

dc_polygons <- SpatialPolygons(list(Polygons(list(Polygon(polygon_coords)), 
                                             ID="DC")))

Polygon

plot(dc_polygons, main="Washington D.C.", col="lightblue", border="blue", 
     xlim = c(-77.3865, -76.8094), ylim = c(38.76835, 39.11431))
points(dmv_air$longitude, dmv_air$latitude, pch=20, col=heat.colors(3)[rank(dmv_air$aqi)])
text(dmv_air$longitude, dmv_air$latitude, labels=dmv_air$name, pos=4, cex=0.7, offset = 0.5)
lines(line_coords, col="blue", lwd=2)

Shapefile and GeoJSON

Shapefile

library(sf)
dc_wards <- st_read("./data/Wards_from_2022/Wards_from_2022.shp")

Reading layer `Wards_from_2022' from data source 
  `/Users/baole/Desktop/spatial-workshop/data/Wards_from_2022/Wards_from_2022.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 8 features and 325 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.1198 ymin: 38.79164 xmax: -76.90915 ymax: 38.99597
Geodetic CRS:  WGS 84

Shapefile and GeoJSON

GEOJSON

airbnb_nb <- st_read("./data/neighbourhoods.geojson")

Reading layer `neighbourhoods' from data source 
  `/Users/baole/Desktop/spatial-workshop/data/neighbourhoods.geojson' 
  using driver `GeoJSON'
Simple feature collection with 39 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.11934 ymin: 38.81801 xmax: -76.90915 ymax: 38.99597
Geodetic CRS:  WGS 84

Projection

Three-dimensional angular system to a two-dimensional planar system.
A coordinate reference system (CRS) defines, with the help of coordinates, how the two-dimensional, projected map is related to real locations on the earth.
- A planar CRS comprises a projection, a datum, and a set of parameters that specify aspects like the map’s center.
The most commonly used datum globally is WGS84 (World Geodesic System 1984), closely related to NAD83 (North American Datum of 1983)
- China: CGCS2000
Most CRSs are assigned an “EPSG code” (European Petroleum Survey Group), which provides a unique identifier to a CRS.
- EPSG:4326, also known as the WGS84 projection (because it’s based on WGS84’s ellipsoid), is a coordinate system used in Google Earth and GSP systems. EPSG:4269 is based on NAD83. EPSG:4490 is based on CGCS2000.

在空间数据处理中，投影（Projection）是一个非常关键的概念。我们通常需要将地球的三维球体表面转换为二维平面，以便在地图上展示空间数据。这一转换过程涉及复杂的数学计算，但其核心是如何精确地将真实地理位置表示在平面上。

什么是投影？投影是将三维地球表面转换为二维地图的一种方法。因为地球是一个椭球体，而地图是平面的，所以投影不可避免地会产生某种程度的失真。投影过程会涉及坐标参考系统（CRS），这是定义投影如何将空间数据映射到地球真实位置的核心。

坐标参考系统 (CRS)：CRS通过坐标定义了二维投影地图与地球上实际位置的关系。CRS通常包括三个要素：投影类型、基准面（Datum）和参数设置（如地图的中心位置）。基准面 (Datum)：基准面是投影的基准点，不同的基准面会影响数据在地图上的精确位置。全球最常用的基准面是WGS84，即世界大地测量系统1984。在北美地区，NAD83（北美基准面1983）也非常常用。 2. 常见的EPSG代码每种CRS都被分配了一个独特的“EPSG代码”。EPSG代码由欧洲石油调查集团（European Petroleum Survey Group）定义，用来唯一标识一种坐标参考系统。

EPSG:4326：这是基于WGS84的坐标系统，也是我们在Google Earth和GPS系统中最常见的坐标系。EPSG:4326中的数据通常使用经纬度（Latitude和Longitude）表示位置。 EPSG:4269：这是基于NAD83的坐标系统，广泛应用于北美地区。对于北美的空间数据应用，使用NAD83能提供更高的地理精度。通过使用标准化的EPSG代码，我们可以在不同的地理信息系统中保持一致的坐标系统，从而确保空间数据的准确性和一致性。

中国常用的投影在中国，常用的投影系统包括CGCS2000（2000年国家大地坐标系）和高斯-克吕格投影（Gauss-Krüger Projection）。

CGCS2000：这是中国的国家大地坐标系，与WGS84类似，但进行了适配和优化，使之更适合中国的地理环境。高斯-克吕格投影：这种投影采用分带系统，即将中国划分为多个投影带，每个带内的数据会使用不同的坐标系。这种方法能够减少投影失真，特别适合大尺度的空间数据，如城市规划和工程测量。

Projection and CRS

head(dc_wards)

Simple feature collection with 6 features and 325 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.08172 ymin: 38.79164 xmax: -76.90915 ymax: 38.9573
Geodetic CRS:  WGS 84
  WARD   NAME          REP_NAME
1    8 Ward 8 Trayon White, Sr.
2    6 Ward 6     Charles Allen
3    7 Ward 7      Vincent Gray
4    2 Ward 2      Brooke Pinto
5    1 Ward 1    Brianne Nadeau
6    5 Ward 5    Zachary Parker
                                                             WEB_URL
1     https://www.dccouncil.us/council/councilmember-trayon-white-sr
2               https://www.dccouncil.us/council/councilmember-allen
3                          https://dccouncil.us/council/vincent-gray
4 https://www.dccouncil.us/council/ward-2-councilmember-brooke-pinto
5                        https://dccouncil.us/council/brianne-nadeau
6                       https://dccouncil.us/council/kenyan-mcduffie
       REP_PHONE             REP_EMAIL
1 (202) 724-8045   twhite@dccouncil.us
2 (202) 724-8072   callen@dccouncil.us
3 (202) 724-8068    vgray@dccouncil.us
4 (202) 724-8058   bpinto@dccouncil.us
5 (202) 724-8181  bnadeau@dccouncil.us
6 (202) 724-8028 zparker@dccouncil.gov
                                  REP_OFFICE WARD_ID  LABEL STUSAB SUMLEV
1 1350 Pennsylvania Ave, Suite 400, NW 20004       8 Ward 8     DC    610
2 1350 Pennsylvania Ave, Suite 110, NW 20004       6 Ward 6     DC    610
3 1350 Pennsylvania Ave, Suite 406, NW 20004       7 Ward 7     DC    610
4 1350 Pennsylvania Ave, Suite 106, NW 20004       2 Ward 2     DC    610
5 1350 Pennsylvania Ave, Suite 108, NW 20004       1 Ward 1     DC    610
6 1350 Pennsylvania Ave, Suite 102, NW 20004       5 Ward 5     DC    610
           GEOID GEOCODE STATE POP100 HU100 P0010001 P0010002 P0010003 P0010004
1 610U600US11008   11008    11  85246 39164    85246    81405     8096    70959
2 610U600US11006   11006    11  84266 52768    84266    77187    50624    19562
3 610U600US11007   11007    11  85685 38968    85685    81893     7038    71628
4 610U600US11002   11002    11  89485 53217    89485    81571    57597    10196
5 610U600US11001   11001    11  85285 45694    85285    76084    42307    18741
6 610U600US11005   11005    11  89617 41794    89617    82833    22248    51307
  P0010005 P0010006 P0010007 P0010008 P0010009 P0010010 P0010011 P0010012
1      305      751       34     1260     3841     3341      988      148
2      216     5028       54     1703     7079     6449      851      388
3      343      591       40     2253     3792     3299      935      140
4      317     9853       64     3544     7914     7407      726      359
5      696     5273       88     8979     9201     8497      956      330
6      495     2731       87     5965     6784     6057     1189      241
  P0010013 P0010014 P0010015 P0010016 P0010017 P0010018 P0010019 P0010020
1      290       22      608      513      156       26      483        4
2     1528       47     2847      228      157       11      274       13
3      219       17      654      591      150       13      486        5
4     1666       34     4056      120       95       20      193       10
5     1294       38     4558      130      197        3      638       10
6      838       22     2146      396      260       12      736        7
  P0010021 P0010022 P0010023 P0010024 P0010025 P0010026 P0010027 P0010028
1        0       46       30       24        3      427      217       32
2        0       31       31       38        5      518      125       73
3        2       32       29       19        7      444      233       42
4        2       61       21       42        2      431       97       51
5        8      239       21       61       14      609      155       50
6        2      133       23       38       14      630      242       40
  P0010029 P0010030 P0010031 P0010032 P0010033 P0010034 P0010035 P0010036
1        5       66        7        1       33       18       20        1
2        4      106       11        0      106       20       33        0
3        1       53        5        0       24       10       21        9
4        2       81        9        0       96       23       34        4
5        5      114       24        0      145       16       49        0
6        1      123        8        1       81       30       42        1
  P0010037 P0010038 P0010039 P0010040 P0010041 P0010042 P0010043 P0010044
1        9        3        5        3        3        0        0        0
2        8        3       17        5        5        1        0        0
3       10        3       15        6        7        2        1        0
4        5        3        9        3        8        4        0        0
5        7        3       19        2        9        6        0        0
6        8        4       25       10        3        9        1        0
  P0010045 P0010046 P0010047 P0010048 P0010049 P0010050 P0010051 P0010052
1        2        2       66       15        3       34        2        6
2        0        1      101       18        0       57        7       14
3        1        1       43       12        0       23        0        3
4        0        2       74       25        0       39        1        2
5        1        4       87       10        1       52        2       10
6        0        1       92       17        0       55        1        7
  P0010053 P0010054 P0010055 P0010056 P0010057 P0010058 P0010059 P0010060
1        0        0        4        0        2        0        0        0
2        0        0        1        0        2        2        0        0
3        0        1        2        0        1        0        0        0
4        0        0        2        0        4        1        0        0
5        0        0        5        0        4        2        1        0
6        0        0        3        3        0        4        1        0
  P0010061 P0010062 P0010063 P0010064 P0010065 P0010066 P0010067 P0010068
1        0        0        7        4        1        0        1        0
2        0        0       10        4        5        0        0        0
3        1        0        6        0        4        0        1        1
4        0        0        1        1        0        0        0        0
5        0        0        8        1        3        0        1        2
6        1        0        5        4        0        0        0        0
  P0010069 P0010070 P0010071 P0020001 P0020002 P0020003 P0020004 P0020005
1        1        0        0    85246     3081    82165    79340     7670
2        1        1        1    84266     6122    78144    73945    49107
3        0        0        0    85685     4078    81607    78876     6624
4        0        1        1    89485     9744    79741    75843    55391
5        1        0        0    85285    17269    68016    64204    40017
6        1        0        0    89617    10419    79198    75285    21217
  P0020006 P0020007 P0020008 P0020009 P0020010 P0020011 P0020012 P0020013
1    70357      213      731       29      340     2825     2475      935
2    19243      121     4977       50      447     4199     3885      800
3    70969      239      568       32      444     2731     2384      907
4    10000      120     9803       54      475     3898     3669      695
5    18330      163     5202       70      422     3812     3506      862
6    50554      218     2687       76      533     3913     3516     1122
  P0020014 P0020015 P0020016 P0020017 P0020018 P0020019 P0020020 P0020021
1      132      272       19      104      482      152       23      297
2      338     1492       40      630      206      153       11      156
3       93      205       11      131      571      143       13      262
4      297     1633       27      662       99       87       20       85
5      238     1279       30      548      118      189        3      172
6      187      802       14      338      381      251       11      349
  P0020022 P0020023 P0020024 P0020025 P0020026 P0020027 P0020028 P0020029
1        4        0       12       29       14        0      315      194
2       12        0        2       22       23        0      268       97
3        3        2        6       29        8        0      326      202
4        9        2        2       15       36        0      200       72
5        8        7        2       19       29        2      275      103
6        7        1        9       22       21        1      368      197
  P0020030 P0020031 P0020032 P0020033 P0020034 P0020035 P0020036 P0020037
1       30        5       37        3        0        7       15        6
2       62        3       38        8        0        5       18       11
3       36        1       35        5        0        4       10        4
4       45        2       23        6        0        3       19        7
5       48        3       48       19        0        4       15        2
6       33        1       53        7        0        1       26       13
  P0020038 P0020039 P0020040 P0020041 P0020042 P0020043 P0020044 P0020045
1        0        8        2        2        3        2        0        0
2        0        7        3        5        5        5        1        0
3        0        8        2        7        6        5        0        1
4        3        4        2        6        3        4        1        0
5        0        6        3       11        2        8        1        0
6        1        6        4        8       10        0        7        1
  P0020046 P0020047 P0020048 P0020049 P0020050 P0020051 P0020052 P0020053
1        0        0        1       32       13        3       10        1
2        0        0        0       40       16        0        9        6
3        0        0        0       19       10        0        7        0
4        0        0        0       28       20        0        2        1
5        0        0        2       27        5        0       10        2
6        0        0        0       27       13        0        8        0
  P0020054 P0020055 P0020056 P0020057 P0020058 P0020059 P0020060 P0020061
1        5        0        0        0        0        0        0        0
2        8        0        0        0        0        0        1        0
3        1        0        1        0        0        0        0        0
4        2        0        0        0        0        2        1        0
5        7        0        0        0        0        1        2        0
6        5        0        0        0        0        0        0        1
  P0020062 P0020063 P0020064 P0020065 P0020066 P0020067 P0020068 P0020069
1        0        0        0        3        2        0        0        1
2        0        0        0        6        4        1        0        0
3        0        0        0        2        0        0        0        1
4        0        0        0        1        1        0        0        0
5        0        0        0        4        1        0        0        1
6        0        0        0        2        2        0        0        0
  P0020070 P0020071 P0020072 P0020073 P0030001 P0030002 P0030003 P0030004
1        0        0        0        0    62507    59754     6972    50933
2        0        1        0        0    74262    68775    45847    16404
3        1        0        0        0    66683    63965     5876    55656
4        0        0        0        0    83628    76770    54915     9062
5        2        0        0        0    75254    67609    39386    15702
6        0        0        0        0    74096    69224    19777    41961
  P0030005 P0030006 P0030007 P0030008 P0030009 P0030010 P0030011 P0030012
1      237      646       25      941     2753     2413      675      111
2      177     4765       50     1532     5487     4988      610      332
3      265      504       33     1631     2718     2375      667       90
4      272     9442       53     3026     6858     6476      630      322
5      543     4954       82     6942     7645     7102      757      291
6      414     2490       75     4507     4872     4360      798      193
  P0030013 P0030014 P0030015 P0030016 P0030017 P0030018 P0030019 P0030020
1      177       19      464      378      104       19      375        4
2      966       42     2422      181      116        7      221       12
3      114       12      498      469       92       11      352        3
4     1372       31     3677       97       61       16      159       10
5     1001       31     3905       94      167        3      554        9
6      480       14     1605      327      181       10      582        6
  P0030021 P0030022 P0030023 P0030024 P0030025 P0030026 P0030027 P0030028
1        0       41       25       19        2      293      161       14
2        0       24       24       30        1      404      101       48
3        1       21       24       14        7      306      159       33
4        1       48       20       30        2      317       73       32
5        6      203       19       49       13      469      129       35
6        2      100       22       28       12      449      181       23
  P0030029 P0030030 P0030031 P0030032 P0030033 P0030034 P0030035 P0030036
1        5       45        4        1       24       10        9        0
2        4       91        4        0       96       11       18        0
3        1       36        2        0       15        6       14        8
4        1       58        7        0       79       16       20        4
5        4       87       15        0      134        8       29        0
6        0       95        3        1       59       23       25        1
  P0030037 P0030038 P0030039 P0030040 P0030041 P0030042 P0030043 P0030044
1        8        3        3        2        2        0        0        0
2        5        3       12        4        5        1        0        0
3        7        2       12        2        7        2        0        0
4        5        3        7        2        5        3        0        0
5        4        2        8        1        4        6        0        0
6        5        2       20        7        0        2        1        0
  P0030045 P0030046 P0030047 P0030048 P0030049 P0030050 P0030051 P0030052
1        0        2       43       10        3       24        1        2
2        0        1       88       15        0       52        5       12
3        0        0       32        7        0       20        0        3
4        0        2       64       24        0       33        1        0
5        0        3       68        4        1       44        2        8
6        0        1       60       11        0       40        0        3
  P0030053 P0030054 P0030055 P0030056 P0030057 P0030058 P0030059 P0030060
1        0        0        1        0        2        0        0        0
2        0        0        1        0        2        1        0        0
3        0        0        1        0        1        0        0        0
4        0        0        2        0        3        1        0        0
5        0        0        3        0        4        2        0        0
6        0        0        2        3        0        1        0        0
  P0030061 P0030062 P0030063 P0030064 P0030065 P0030066 P0030067 P0030068
1        0        0        4        3        1        0        0        0
2        0        0        7        3        4        0        0        0
3        0        0        5        0        4        0        0        1
4        0        0        1        1        0        0        0        0
5        0        0        6        1        3        0        0        2
6        0        0        3        3        0        0        0        0
  P0030069 P0030070 P0030071 P0040001 P0040002 P0040003 P0040004 P0040005
1        0        0        0    62507     2195    60312    58274     6649
2        0        0        0    74262     5257    69005    65905    44499
3        0        0        0    66683     2938    63745    61796     5529
4        0        0        0    83628     8618    75010    71714    52844
5        0        0        0    75254    13901    61353    58271    37464
6        0        0        0    74096     7776    66320    63538    18976
  P0040006 P0040007 P0040008 P0040009 P0040010 P0040011 P0040012 P0040013
1    50546      178      628       23      250     2038     1779      651
2    16156      102     4722       48      378     3100     2854      581
3    55244      197      484       25      317     1949     1718      648
4     8908      111     9397       46      408     3296     3121      605
5    15389      129     4890       65      334     3082     2860      690
6    41437      180     2453       69      423     2782     2492      763
  P0040014 P0040015 P0040016 P0040017 P0040018 P0040019 P0040020 P0040021
1       99      168       17       83      358      100       16      241
2      303      947       36      514      169      114        7      131
3       65      106       10       98      454       86       11      201
4      271     1344       26      585       83       59       16       80
5      218      993       24      473       86      163        3      150
6      151      462        7      261      317      175        9      292
  P0040022 P0040023 P0040024 P0040025 P0040026 P0040027 P0040028 P0040029
1        4        0        7       25       10        0      233      152
2       11        0        2       18       21        0      207       82
3        1        1        5       24        8        0      218      136
4        9        1        1       15       26        0      150       56
5        7        5        2       18       26        2      203       88
6        6        1        9       21       17        1      273      157
  P0040030 P0040031 P0040032 P0040033 P0040034 P0040035 P0040036 P0040037
1       14        5       26        3        0        6        7        4
2       40        3       34        4        0        5       10        9
3       27        1       21        2        0        2        6        2
4       26        1       17        5        0        3       15        7
5       34        2       37       12        0        3        8        2
6       17        0       45        3        0        1       20        9
  P0040038 P0040039 P0040040 P0040041 P0040042 P0040043 P0040044 P0040045
1        0        7        2        2        2        2        0        0
2        0        4        3        3        4        5        1        0
3        0        7        1        6        2        5        0        0
4        3        4        2        5        2        4        0        0
5        0        3        2        4        1        4        1        0
6        1        3        2        7        7        0        0        1
  P0040046 P0040047 P0040048 P0040049 P0040050 P0040051 P0040052 P0040053
1        0        0        1       25        9        3       10        1
2        0        0        0       35       13        0        8        5
3        0        0        0       12        6        0        5        0
4        0        0        0       24       20        0        1        1
5        0        0        2       16        1        0        5        2
6        0        0        0       16        7        0        6        0
  P0040054 P0040055 P0040056 P0040057 P0040058 P0040059 P0040060 P0040061
1        2        0        0        0        0        0        0        0
2        8        0        0        0        0        0        1        0
3        1        0        0        0        0        0        0        0
4        0        0        0        0        0        1        1        0
5        5        0        0        0        0        1        2        0
6        3        0        0        0        0        0        0        0
  P0040062 P0040063 P0040064 P0040065 P0040066 P0040067 P0040068 P0040069
1        0        0        0        1        1        0        0        0
2        0        0        0        4        3        1        0        0
3        0        0        0        1        0        0        0        0
4        0        0        0        1        1        0        0        0
5        0        0        0        3        1        0        0        0
6        0        0        0        1        1        0        0        0
  P0040070 P0040071 P0040072 P0040073 H0010001 H0010002 H0010003 P0050001
1        0        0        0        0    39164    34608     4556     3362
2        0        0        0        0    52768    45285     7483     1991
3        1        0        0        0    38968    34958     4010     4528
4        0        0        0        0    53217    46543     6674    13218
5        2        0        0        0    45694    41035     4659     5066
6        0        0        0        0    41794    37211     4583     6759
  P0050002 P0050003 P0050004 P0050005 P0050006 P0050007 P0050008 P0050009
1     1054       23       32      750      249     2308        0      563
2      194       31       29      134        0     1797      655      255
3     2557     2197       34      289       37     1971        0        0
4      205       27        0      178        0    13013    11470        0
5      305        0      173      132        0     4761     4009        0
6     1007        0       40      967        0     5752     3053        0
  P0050010 OBJECTID                               GLOBALID CREATED_US
1     1745        1 {E31550AE-6FAE-4B74-909F-52B283BFAF68}       <NA>
2      887        2 {765C4F49-9292-4BDB-AA24-39F4EE43359F}       <NA>
3     1971        3 {73F07042-7D7F-452B-9BB3-0F87B0EC5418}       <NA>
4     1543        4 {7F8C2A51-427C-45FC-91EB-9693656AED9C}       <NA>
5      752        5 {C3C6E2E7-E68D-49B2-970C-D60675EA7B4B}       <NA>
6     2699        6 {6C10DD95-DE70-4F26-94BA-79F378FA74E0}       <NA>
  CREATED_DA LAST_EDITE LAST_EDI_1 SHAPEAREA SHAPELEN
1       <NA>       <NA>       <NA>         0        0
2       <NA>       JLAY 2023-12-07         0        0
3       <NA>       <NA>       <NA>         0        0
4       <NA>       <NA>       <NA>         0        0
5       <NA>       JLAY 2023-12-07         0        0
6       <NA>       JLAY 2023-12-07         0        0
                        geometry
1 POLYGON ((-76.99392 38.8777...
2 POLYGON ((-77.00908 38.8716...
3 POLYGON ((-76.94186 38.9185...
4 POLYGON ((-77.0327 38.88262...
5 POLYGON ((-77.03523 38.9374...
6 POLYGON ((-76.99144 38.9573...

Projection and CRS

st_crs(airbnb_nb)

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["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Spatial Weights Matrix

A vast number of spatial analysis methods are based on a spatial matrix \mathbf{W}_{ij} of size n \times n with i an observation and j the “neighbours” of that observation. \mathbf{W}_{ij} represents the degree of spatial relationship between i and j.
Classically one can define :
- Neighbouring matrix (\mathbf{W}_{ij} = 1 if j is a neighbour of i, 0 otherwise)
- Distance matrix (\mathbf{W}_{ij} = the distance between i and j, modified by a function like \frac{1}{distance} or \frac{1}{distance^2})
- Interaction matrix (\mathbf{W}_{ij} = the degree of interaction between i and j, the measure of the interaction depends on the subject of the analysis)

在空间分析中，空间权重矩阵（Spatial Weights Matrix）是一个核心概念。很多空间分析方法都依赖于一个n \times n的权重矩阵\mathbf{W}{ij}，其中i和j分别代表观察点及其邻居。\mathbf{W}{ij}描述了观察点i与邻居j之间的空间关系强度。这种关系的表示方式可以根据分析需求的不同而有所差异。

什么是空间权重矩阵？空间权重矩阵\mathbf{W}_{ij}定义了空间单元之间的关联关系。在空间分析中，常用的矩阵类型包括邻接矩阵、距离矩阵和交互矩阵。通过这些不同的矩阵，我们可以灵活地定义空间单元之间的关系：

邻接矩阵：最简单的空间权重矩阵类型。通常，如果j是i的邻居，则\mathbf{W}{ij} = 1；如果不是，则\mathbf{W}{ij} = 0。这种矩阵用于表示地理上相邻的空间单元，是分析空间相邻效应的基础。

距离矩阵：在距离矩阵中，\mathbf{W}_{ij}表示i和j之间的距离。距离可以通过函数进行调整，比如$ 或 $，以反映距离对空间影响的衰减效果。这种方法特别适合用于分析随着距离变化而减弱的空间效应，例如空气污染的扩散。

交互矩阵：交互矩阵则是根据具体分析主题来定义的。\mathbf{W}_{ij}表示i和j之间的交互强度，比如人口流动、贸易量等。这种矩阵非常灵活，可以根据研究需求来量化特定类型的空间关系。

空间权重矩阵的应用空间权重矩阵在空间统计和地理分析中有广泛应用。通过构建权重矩阵，我们能够捕捉空间自相关性，即地理上相邻的空间单元在行为和特征上的相似性。空间权重矩阵的类型选择将直接影响分析结果。例如：

邻接矩阵在房价、犯罪率等相邻效应研究中常用，因为它能有效描述相邻区域之间的相似性。距离矩阵常用于环境科学和公共卫生研究，通过距离衰减效应分析污染、疾病传播等现象。交互矩阵则广泛应用于社会经济研究，比如流动人口的空间影响分析，通过交互强度来反映单元之间的联系程度。总结空间权重矩阵\mathbf{W}_{ij}是空间分析的基础工具。它灵活地描述了空间单元之间的关系强度，不同的矩阵类型能够满足多种研究需求。理解空间权重矩阵的结构和定义方式，将有助于我们更深入地进行空间数据分析和地理关系建模。

Spatial Weights Matrix: Neighborhood

library(spdep)
library(sf)
columbus <- st_read(system.file("shapes/columbus.shp",
               package = "spData"), quiet = TRUE)
nb <- poly2nb(columbus, queen = TRUE)
head(nb)

[[1]]
[1] 2 3

[[2]]
[1] 1 3 4

[[3]]
[1] 1 2 4 5

[[4]]
[1] 2 3 5 8

[[5]]
[1]  3  4  6  8  9 11 15 16

[[6]]
[1] 5 9

Spatial Neighborhood

车邻接 (Rook Contiguity) 车邻接规则类似于国际象棋中的车移动方式，要求两个单元必须共享边界（不包括顶点）才能视为邻居。

特点：车邻接相对更严格，仅允许共享边界的单元相邻，排除了仅通过顶点相接的关系。应用场景：车邻接更适合用于分析共享边界的区域，比如分区管理、生态研究中的栖息地边界等。代码示例：在R中，通过将queen参数设为FALSE，我们可以生成车邻接矩阵。

主教邻接 (Bishop Contiguity) 主教邻接定义的是仅通过顶点相邻的关系，类似于国际象棋中的主教移动方式。

特点：这种邻接规则只包括顶点相邻的单元，而排除了边界相邻的关系，因此较为特殊。应用场景：主教邻接在空间分析中应用较少，但在特定的地理空间模型中可能有用，例如捕捉特定角落邻接的区域。代码实现：主教邻接可以通过手动筛选顶点邻接的单元来实现，R中的常见包不直接支持此规则，需要自定义代码。

皇后邻接 (Queen Contiguity) 皇后邻接规则类似于国际象棋中的皇后移动方式，定义两个单元相邻的条件是它们共享边界或顶点。

特点：这种邻接规则较为宽松，允许单元通过顶点相邻，因此可以捕捉更多的空间联系。应用场景：皇后邻接适用于那些邻居之间互动较强的情况，比如城市中的社区或街区分析。代码示例：在R中，我们可以通过poly2nb(…, queen = TRUE)来构建皇后邻接矩阵。

六边形邻接 (Hexagon Neighborhood) 六边形邻接是一种专门用于六边形网格的邻接规则，因为六边形的几何结构特点，每个六边形都有六个邻居。

特点：六边形网格在空间分析中具有均匀的邻接结构，避免了矩形网格中边界长度不一致的问题，因此六边形网格常用于生态学和环境建模。应用场景：六边形邻接适用于连续空间的建模，特别是在环境监测、生态模拟和地理数据采样中被广泛使用。代码示例：在R中可以使用spdep包的hexagonal功能创建六边形网格。

Neighborhood

plot(st_geometry(columbus), border = "lightgray")
plot(nb, st_geometry(columbus), add = TRUE)

Neighborhood

id <- 25 # select an area
columbus$neighbors <- "other"
columbus$neighbors[id] <- "area"
columbus$neighbors[nb[[id]]] <- "neighbors"
cols <- c("gray30", "gray", "white")
plot(st_geometry(columbus), col=cols[as.factor(columbus$neighbors)], main="Spatial Neighborhoods")
legend("topright", legend=c("Selected Area", "Neighbors", "Other"), fill=cols, bty="n")

Neighborhood

Neighborhood Matrix

k Nearest Neighbors

E.g.: neighbors based on 3 nearest neighbors

k Nearest Neighbors

cols_ctrs <- st_centroid(columbus)
knb <- knn2nb(knearneigh(cols_ctrs, k = 5)) # k=3 number nearest neighbors
plot(st_geometry(columbus), border = "lightgray")
plot(knb, st_geometry(columbus), add = TRUE)

k近邻是一种基于距离的邻接方式。对于每个地理单元i，我们根据地理距离从近到远为其选择k个邻居单元。这样，无论空间单元的形状或位置如何，每个单元都能拥有固定数量的邻居。

灵活性：k-NN不依赖于地理单元的形状，因此在空间布局不规则的情况下非常有效，比如分析散点分布的数据。参数k的选择：k的大小通常根据具体的分析需求来确定。较小的k值可以捕捉更局部的关系，而较大的k值则适合捕捉更广泛的空间关系。

k近邻方法不仅在地理空间分析中使用广泛，在其他领域也非常流行，尤其是在机器学习和数据挖掘中。以下是一些典型的应用示例：

分类问题：在机器学习中，k-NN算法是一种常用的分类算法。对于一个未知类别的样本，k-NN会根据它与已知样本的距离，将它归入与它最接近的k个邻居中出现频率最高的类别。这种方法在图像识别、文本分类等领域都有广泛应用。

推荐系统：在推荐系统中，k-NN用于查找用户或产品的“相似邻居”。例如，Netflix可以利用k-NN为用户推荐相似的电影或电视剧，基于其他具有相似观看记录的用户偏好。

k-NN在模式识别中非常常用，尤其是在医学图像处理和基因序列分析中。比如，k-NN可以根据相似的特征模式来识别不同的细胞类型或基因序列。

Spatial Weight Matrix: Distance

Distance matrix

dist_mat <- spDists(as.matrix(st_coordinates(cols_ctrs)), longlat = FALSE)
dist_mat[1:6,1:6]

          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
[1,] 0.0000000 0.5987183 0.5796856 0.7480069 1.0878333 1.2779137
[2,] 0.5987183 0.0000000 0.7118774 0.3397537 0.9983315 1.5134864
[3,] 0.5796856 0.7118774 0.0000000 0.5609564 0.5233702 0.8053411
[4,] 0.7480069 0.3397537 0.5609564 0.0000000 0.6901508 1.3018772
[5,] 1.0878333 0.9983315 0.5233702 0.6901508 0.0000000 0.7533395
[6,] 1.2779137 1.5134864 0.8053411 1.3018772 0.7533395 0.0000000

Inverse distance weights (IDW) matrix

inv_dist_mat <- 1 / (dist_mat)
diag(inv_dist_mat) <- 0
inv_dist_mat[1:6,1:6]

          [,1]      [,2]     [,3]      [,4]      [,5]      [,6]
[1,] 0.0000000 1.6702347 1.725073 1.3368861 0.9192585 0.7825254
[2,] 1.6702347 0.0000000 1.404736 2.9433086 1.0016712 0.6607261
[3,] 1.7250730 1.4047363 0.000000 1.7826697 1.9106935 1.2417099
[4,] 1.3368861 2.9433086 1.782670 0.0000000 1.4489588 0.7681216
[5,] 0.9192585 1.0016712 1.910694 1.4489588 0.0000000 1.3274227
[6,] 0.7825254 0.6607261 1.241710 0.7681216 1.3274227 0.0000000

Inverse Distance Weights

id <- 25
columbus$weights <- inv_dist_mat[id, ]

# Plot the map, coloring regions based on the selected area's weights
plot(columbus["weights"], col = gray(1 - columbus$weights/max(columbus$weights)), main = "Inverse Distance Weights")

Autocorrelation

Spatial autocorrelation is used to describe the extent to which a variable is correlated with itself through space.
- Positive spatial autocorrelation occurs when observations with similar values are closer together (i.e., clustered).
Global vs. local spatial autocorrelation
Statistical tests
- Moran’s I (Moran 1950): correlation coefficient.
- Geary’s C (Geary 1954): distance-based and more sensitive to local spatial autocorrelation.

在空间分析中，空间自相关（Spatial Autocorrelation）是描述变量在地理空间上自相关程度的关键概念。空间自相关帮助我们理解某一变量的相似值是否倾向于在空间上聚集，或是随机分布。接下来，我们深入探讨空间自相关的含义、类型以及常用的统计检验方法。

什么是空间自相关？空间自相关用于描述一个变量在空间上与自身的关联程度。当相似的观测值聚集在一起时，我们称之为正空间自相关。例如，在城市中，高收入家庭往往集中在某些特定区域，而低收入家庭集中在其他区域，这就是典型的正空间自相关现象。

正空间自相关：当相似的值在空间上靠近时，空间自相关为正，表现为聚集性模式。负空间自相关：当相似值分布在更远的空间上时，空间自相关为负，表现为离散模式。 2. 全局与局部空间自相关空间自相关可以分为全局和局部空间自相关：

全局空间自相关：全局自相关测量变量在整个空间范围内的自相关程度，用于判断是否存在整体性的聚集或离散模式。全局自相关能够揭示广泛的空间趋势。局部空间自相关：局部自相关则聚焦于局部区域，识别小范围内的空间模式。例如，局部自相关可以帮助我们识别特定的“热点”区域或“冷点”区域，从而深入了解空间结构的细节。

常用的统计检验方法在空间自相关分析中，Moran’s I和Geary’s C是两种常用的统计检验方法。

Moran’s I：由Moran在1950年提出，是一种基于相关系数的统计方法，用于测量全局空间自相关。Moran’s I的值介于-1和1之间，正值表示正空间自相关，负值表示负空间自相关。Moran’s I在研究全局趋势和广泛的空间分布模式时非常有效。

Geary’s C：由Geary在1954年提出，是一种基于距离的统计方法。Geary’s C对局部空间自相关更为敏感，能够捕捉到小范围的空间差异。与Moran’s I不同，Geary’s C的值越接近0，表示正空间自相关；越接近2，则表示负空间自相关。

Autocorrelation: Moran’s I

Global Moran’s I

I = \frac{n \sum\limits_i \sum\limits_j w_{ij}(Y_i - \bar{Y})(Y_j - \bar{Y})}{\left(\sum\limits_{i \neq j} w_{ij}\right) \sum\limits_i (Y_i - \bar{Y})^2}

Expected value of Moran’s I: E(I) = \frac{-1}{n-1}
Test statistic: z = \frac{I - E(I)}{\sqrt{Var(I)}}

Moran’s I

To what extent the Moran’s I value is different from random

Moran’s I

Local Moran’s I

I_i = \frac{n(Y_i - \bar{Y})}{\sum_j (Y_j - \bar{Y})^2} \sum_j w_{ij}(Y_j - \bar{Y})

Note that the global Moran’s I is proportional to the sum of the local Moran’s I obtained for all locations:

I = \frac{1}{\sum\limits_{i \neq j} w_{ij}} \sum\limits_{i} I_i

Script for Slide: Moran’s I - Local Moran’s I 除了全局 Moran’s I 之外，我们还有局部 Moran’s I，用于识别特定位置的空间自相关。局部 Moran’s I 能够揭示每个观测单元在空间聚集或离散中的作用，帮助我们识别数据中的“热点”（高值聚集）或“冷点”（低值聚集）。

局部 Moran’s I 的公式在这个公式中： I_i 表示位置 i 的局部 Moran’s I 值。 Y_i 和 Y_j 分别是位置 i 和 j 的属性值。 \bar{Y} 是所有位置的平均值。 w_{ij} 是空间权重矩阵的元素，用于表示 i 和 j 的空间关系强度。局部 Moran’s I 通过对比每个位置的值和邻近位置的平均值来判断该位置的空间自相关性。例如，高值与高值相邻或低值与低值相邻会得到正的局部 Moran’s I，而高值与低值相邻会得到负的局部 Moran’s I。
全局 Moran’s I 与局部 Moran’s I 的关系需要注意的是，全局 Moran’s I 与所有位置的局部 Moran’s I 之间存在比例关系：这个公式表明，全局 Moran’s I 可以看作所有局部 Moran’s I 的加权平均。因此，通过计算每个位置的局部 Moran’s I，我们可以构建出全局 Moran’s I，从而了解整体和局部的空间自相关性。
局部 Moran’s I 的应用局部 Moran’s I 常用于识别空间数据中的“热点”和“冷点”：

热点（High-High Cluster）：某个位置的值高，并且与高值的邻居相邻时，会形成热点区域，代表一个高度聚集的区域。冷点（Low-Low Cluster）：某个位置的值低，并且与低值的邻居相邻时，会形成冷点区域。高-低离散或低-高离散：当高值与低值相邻时，局部 Moran’s I 为负，表示一种离散模式。通过局部 Moran’s I，我们可以发现地理数据中的聚集模式和异常区域，为政策决策和空间规划提供数据支持。

Autocorrelation: Boston House Price

boston <- st_read(system.file("shapes/boston_tracts.shp",
                           package = "spData"), quiet = TRUE)
head(boston[,1:7])

Simple feature collection with 6 features and 7 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -71.1753 ymin: 42.33 xmax: -71.1238 ymax: 42.3737
Geodetic CRS:  NAD27
  poltract                    TOWN TOWNNO TRACT      LON     LAT MEDV
1     0001 Boston Allston-Brighton     74     1 -71.0830 42.2172 17.8
2     0002 Boston Allston-Brighton     74     2 -71.0950 42.2120 21.7
3     0003 Boston Allston-Brighton     74     3 -71.1007 42.2100 22.7
4     0004 Boston Allston-Brighton     74     4 -71.0930 42.2070 22.6
5     0005 Boston Allston-Brighton     74     5 -71.0905 42.2033 25.0
6     0006 Boston Allston-Brighton     74     6 -71.0865 42.2100 19.9
                        geometry
1 POLYGON ((-71.1238 42.3689,...
2 POLYGON ((-71.1546 42.3573,...
3 POLYGON ((-71.1685 42.3601,...
4 POLYGON ((-71.15391 42.3461...
5 POLYGON ((-71.1479 42.337, ...
6 POLYGON ((-71.1382 42.3535,...

Boston House Price

plot(boston["MEDV"], main="Boston median prices of owner-occupied housing ($1000 USD)")

Boston House Price

mapview::mapview(boston, zcol = "MEDV")

Autocorrelation of House Price

nb <- poly2nb(boston) 
nbw <- nb2listw(nb, style = "W")

# Global Moran's I
gmoran <- moran.test(boston$MEDV, nbw,
                     alternative = "greater")
gmoran


    Moran I test under randomisation

data:  boston$MEDV  
weights: nbw    

Moran I statistic standard deviate = 23.35, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.6266753872     -0.0019801980      0.0007248686

Monte Carlo Method of Moran’s I

A permutation bootstraping approach

gmoranMC <- moran.mc(boston$MEDV, nbw, nsim = 1000)
gmoranMC


    Monte-Carlo simulation of Moran I

data:  boston$MEDV 
weights: nbw  
number of simulations + 1: 1001 

statistic = 0.62668, observed rank = 1001, p-value = 0.000999
alternative hypothesis: greater

我们可以使用蒙特卡洛方法（Monte Carlo Method）。蒙特卡洛方法是一种置换引导法（permutation bootstrapping），它通过大量随机置换来计算Moran’s I的分布，从而获得更稳健的显著性检验结果。接下来，我们将详细介绍蒙特卡洛方法在Moran’s I中的应用原理和流程。

什么是蒙特卡洛方法？蒙特卡洛方法是一种模拟技术，通过随机置换数据来估计统计量的分布。在Moran’s I的显著性检验中，蒙特卡洛方法可以帮助我们更可靠地评估观察到的Moran’s I值是否显著高于随机分布。

置换引导法：我们将空间数据中的变量值随机置换多次（例如1000次），每次重新计算Moran’s I。这样可以建立一个在随机条件下的Moran’s I分布，从而与实际观测值进行比较。显著性评估：通过比较观察值与随机置换产生的Moran’s I分布，我们可以判断观测到的Moran’s I是否显著。 2. 蒙特卡洛方法的流程蒙特卡洛方法应用在Moran’s I显著性检验中的具体流程如下：

随机置换：对变量值进行多次随机置换，确保每次置换后空间单元的值分布不同。计算Moran’s I：每次置换后重新计算Moran’s I，生成多个（例如1000次）随机条件下的Moran’s I值。建立分布：将所有置换后的Moran’s I值构成一个随机分布。显著性检验：将实际观测的Moran’s I值与随机分布进行比较，计算其在分布中的百分位位置。如果观测值位于分布的高端或低端，则说明Moran’s I显著不同于随机分布。 3. 结果解释如果观测到的Moran’s I值显著高于随机置换分布中的大多数值，则我们可以认为存在显著的正空间自相关性。相反，如果观察值显著低于大多数置换值，则表示存在负空间自相关性。通过蒙特卡洛方法，Moran’s I显著性检验的结果更加稳健，特别适合用于复杂或不规则空间数据的分析。

Monte Carlo Method of Moran’s I

hist(gmoranMC$res, breaks = round(seq(min(gmoranMC$res)-0.01,max(gmoranMC$res)+0.01,0.01), 3), 
     main = "Distribution of Results")
abline(v = gmoranMC$statistic, col = "red", lwd=1.5, lty=2)

Clusters Based on Autocorrelation

We can use local Moran’s I to identify cluster.

lmoran <- localmoran(boston$MEDV, nbw, alternative = "two.sided")
head(lmoran)

             Ii          E.Ii       Var.Ii        Z.Ii Pr(z != E(Ii))
1 -0.3457508492 -5.254157e-04 3.275376e-02 -1.90753363   5.645152e-02
2  0.0175875407 -1.626873e-05 2.045711e-03  0.38921049   6.971204e-01
3  0.0123379633 -6.557001e-07 4.089699e-05  1.92939381   5.368199e-02
4 -0.0001654033 -1.059064e-07 1.331742e-05 -0.04529559   9.638717e-01
5  0.3591628595 -1.427815e-04 7.898947e-03  4.04277384   5.282257e-05
6  0.0545610965 -1.625936e-04 1.357382e-02  0.46970410   6.385664e-01

boston$pvals <- lmoran[, 5]

Clusters

moran.plot(as.vector(boston$MEDV), nbw)
text(12.5, 45, "Low-High"); text(35, 45, "High-High");text(12.5, 7.5, "Low-Low"); text(35, 7.5, "High-Low")

Spatial Interpolation

Spatial dependence: bug or feature
Estimating or predicting values at unsampled locations using observations from sampled locations based on spatial dependence.
- Creating continuous surfaces

空间插值（Spatial Interpolation）是一种在地理空间中估算未采样位置数据的方法。利用空间依赖性，我们可以根据已有采样点的数据，预测其他位置的值，从而生成连续的空间数据表面。空间插值在气象、地质、环境监测等领域应用广泛。接下来，我们将详细探讨空间插值的基本概念和应用。

空间依赖性：障碍还是优势？在空间分析中，空间依赖性既可以是一个分析挑战，也可以是一个强大的工具。在空间插值中，空间依赖性被视为一种“特性”，而非“障碍”。这意味着我们假定邻近位置的数据更相似，从而可以利用这种依赖性进行数据预测。

例子：在气温预测中，空间依赖性表明相邻地点的气温更接近，因此可以通过已知地点的气温估算邻近未采样地点的温度。 2. 空间插值的核心目标空间插值的主要目标是根据已有样本位置的数据，估算未采样位置的值。这种方法可以帮助我们生成一个连续的空间表面，使得数据在地理空间上更加完整和易于分析。

连续表面：例如，在空气质量监测中，我们通常只有若干传感器点的数据。通过空间插值，我们可以生成一个覆盖整个区域的空气质量分布图，为公共卫生分析提供依据。预测新位置：空间插值可以用于预测新位置的值，从而弥补采样点不足的缺陷。即使未采样点没有直接观测值，我们也可以通过邻近点的已知值进行合理估算。 3. 应用场景空间插值广泛应用于各种需要预测或生成连续空间数据的领域：

气象学：通过空间插值，气象学家可以预测不同区域的温度、降水量等气象数据，为天气预报提供支持。地质学：在矿产资源评估中，空间插值可以帮助估算未勘探区域的矿物含量，从而指导资源开采。环境监测：在污染监测中，空间插值生成的污染分布图可以帮助决策者识别高污染区域，制定有针对性的政策。

Spatial Interpolation Methods

“Closest observation” (naive estimates)
Thiessen Polygons

在空间插值中，我们有多种方法可以估算未采样位置的数据值。不同的方法适用于不同的场景，今天我们将介绍其中两种基础插值方法：最近观测值（Closest Observation）和Thiessen多边形。

最近观测值 (Closest Observation) 最简单的空间插值方法是“最近观测值”法。这种方法直接使用离目标位置最近的观测值来估算目标位置的值。

原理：最近观测值方法假定空间依赖性较强，因此未采样位置的数据值可以直接用离它最近的已知观测值来代替。优缺点：这种方法简单且易于计算，但它没有考虑周围多个观测值的影响，因此可能导致预测结果在空间上不连续。最近观测值法适用于需要快速、简单估算的情况，但在精度要求较高的分析中，其适用性有限。

Thiessen 多边形 (Thiessen Polygons) Thiessen 多边形是一种基于空间划分的插值方法，它将空间区域划分为多个多边形，每个多边形中所有点的值等于其中心点的观测值。

原理：每个观测点的影响范围被划分成一个多边形，称为Thiessen多边形或Voronoi多边形。多边形内的所有位置都采用该观测点的值，因为这些位置到该观测点的距离最近。图示：右侧的图展示了Thiessen多边形的划分情况。每个多边形的边界是等距离于相邻观测点的分界线。优缺点：Thiessen多边形方法比最近观测值法更精确，因为它考虑了每个观测点的空间范围，但它仍然没有平滑的过渡效果，导致插值结果在多边形边界处出现不连续的跳变。 Thiessen多边形法适用于土地分配、资源评估等场景中，这些场景中空间内每个位置都依赖于离它最近的一个观测点的值。

Spatial Interpolation Methods

Inverse Distance Weighting (IDW)

\hat{y}(s_0) = \frac{\sum_{i=1}^{n} y(s_i) \times (1/d_i^{\beta})}{\sum_{i=1}^{n} (1/d_i^{\beta})} = \sum_{i=1}^{n} y(s_i) \times w_i

Nearest neighbors

\hat{y}(s_0) = \frac{\sum_{i=1}^{k} y(s_i)}{k}

Ensemble approach
- New development: ensemble of multiple data sources (monitor+satellite+simulation); hyper local; etc.

在空间插值中，我们可以使用多种方法来估算未采样位置的值。这些方法各有特点，适用于不同的应用场景。今天我们将讨论三种方法：反距离权重法（Inverse Distance Weighting, IDW）、最近邻平均法（Nearest Neighbors），以及集成方法（Ensemble Approach）。

反距离权重法 (Inverse Distance Weighting, IDW) 反距离权重法是一种基于距离衰减的插值方法，假定距离越近的观测点对目标位置的影响越大。

原理：在IDW方法中，每个观测值y(s_i)的影响力与它离目标位置的距离d_i成反比，通过指数参数\beta来控制距离衰减的程度。较大的\beta值表示更快速的衰减，强调近邻观测点。优点和应用：IDW能够生成相对平滑的空间表面，适合用于温度、降水量等随距离变化逐渐衰减的自然现象。IDW方法平衡了多个邻近观测点的影响，通常在环境和生态领域中广泛应用。

最近邻平均法 (Nearest Neighbors) 最近邻平均法是一种简单的插值方法，基于选取目标位置周围的k个最近邻居的均值来估算目标位置的值。

原理：对于每个未采样点s_0，最近邻平均法选取k个距离最近的观测点的值，并计算它们的均值作为插值结果。此方法假定最近的观测值足以代表目标位置的情况。优点和应用：最近邻平均法简单易用，适用于需要快速估算的场景。它通常被用于小区域数据或需要低计算成本的场景，但它忽略了距离的渐变影响，可能导致空间不连续性。 1. 集成方法 (Ensemble Approach) 集成方法是一种新兴发展的插值方法，它将来自多个数据源的观测结果进行集成，以获得更精准和“超本地化”（hyper-local）的预测。

原理：集成方法结合了不同来源的数据，例如监测站数据、卫星遥感数据和数值模拟数据。通过综合多种信息来源，集成方法可以提供更准确的空间预测。优点和应用：集成方法特别适用于需要高精度和超本地化预测的场景。例如，空气质量预测中可以集成地面监测数据、卫星图像和气象模拟结果，以便提供更精准的污染分布图。这种方法对环境监测、城市规划等领域的高分辨率预测需求非常有帮助。

Spatial Interpolation: DMV Air Quality

library(ggplot2)
us_county <- st_read("./data/cb_2018_us_county_500k/cb_2018_us_county_500k.shp", quiet=T)
dmv_county <- us_county[us_county$STATEFP %in% c("11", "24", "51"),]
dmv_county

Simple feature collection with 158 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -83.67539 ymin: 36.54074 xmax: -75.04894 ymax: 39.72304
Geodetic CRS:  NAD83
First 10 features:
    STATEFP COUNTYFP COUNTYNS       AFFGEOID GEOID                 NAME LSAD
58       11      001 01702382 0500000US11001 11001 District of Columbia   00
145      24      003 01710958 0500000US24003 24003         Anne Arundel   06
296      51      015 01480098 0500000US51015 51015              Augusta   06
297      51      037 01492442 0500000US51037 51037            Charlotte   06
298      51      051 01497376 0500000US51051 51051            Dickenson   06
299      51      053 01690739 0500000US51053 51053            Dinwiddie   06
300      51      069 01480124 0500000US51069 51069            Frederick   06
301      51      073 01480126 0500000US51073 51073           Gloucester   06
302      51      085 01480132 0500000US51085 51085              Hanover   06
303      51      099 01480137 0500000US51099 51099          King George   06
         ALAND    AWATER                       geometry
58   158340391  18687198 MULTIPOLYGON (((-77.11976 3...
145 1074287145 448099595 MULTIPOLYGON (((-76.84036 3...
296 2504708478   9537636 MULTIPOLYGON (((-79.53328 3...
297 1231010119   5809378 MULTIPOLYGON (((-78.90459 3...
298  855860629   8082181 MULTIPOLYGON (((-82.55383 3...
299 1305091846   8954562 MULTIPOLYGON (((-77.90034 3...
300 1070101640   5935153 MULTIPOLYGON (((-78.54022 3...
301  564117808 181754414 MULTIPOLYGON (((-76.7063 37...
302 1211818312  12652671 MULTIPOLYGON (((-77.78648 3...
303  465122891  21198264 MULTIPOLYGON (((-77.34744 3...

DMV Air Quality

dmv_air <- st_read("./data/dmv_air/dmv_air.shp", quiet=T)
dmv_air

Simple feature collection with 35 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -82.1641 ymin: 36.608 xmax: -75.79732 ymax: 39.70595
Geodetic CRS:  WGS 84
First 10 features:
                Stat_Nm Site_Nm                                    Lcl_S_N
1  District Of Columbia      41                              RIVER TERRACE
2  District Of Columbia      43                       MCMILLAN  NCore-PAMS
3  District Of Columbia      51                                  Near Road
4  District Of Columbia      53                  King Greenleaf Rec Center
5              Maryland       2                                 Millington
6              Maryland       2                                  Piney Run
7              Maryland       3 Fair Hill Natural Resource Management Area
8              Maryland       4                                 Horn Point
9              Maryland       6                    Howard County Near Road
10             Maryland       9                                 Hagerstown
    annl_mn                   geometry
1  9.661741 POINT (-76.95807 38.89557)
2  8.073839 POINT (-77.01318 38.92185)
3  8.695211 POINT (-76.95343 38.89477)
4  8.934325 POINT (-77.01282 38.87516)
5  6.922022 POINT (-75.79732 39.30502)
6  6.403240   POINT (-79.012 39.70595)
7  7.459565 POINT (-75.86477 39.70298)
8  6.918579 POINT (-76.14101 38.58752)
9  8.705889 POINT (-76.84611 39.14313)
10 7.068597 POINT (-77.72024 39.56418)

DMV Air Quality

dmv_air <- st_transform(dmv_air, crs = crs(dmv_county))
dmv_air

Simple feature collection with 35 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -82.1641 ymin: 36.608 xmax: -75.79732 ymax: 39.70595
Geodetic CRS:  NAD83
First 10 features:
                Stat_Nm Site_Nm                                    Lcl_S_N
1  District Of Columbia      41                              RIVER TERRACE
2  District Of Columbia      43                       MCMILLAN  NCore-PAMS
3  District Of Columbia      51                                  Near Road
4  District Of Columbia      53                  King Greenleaf Rec Center
5              Maryland       2                                 Millington
6              Maryland       2                                  Piney Run
7              Maryland       3 Fair Hill Natural Resource Management Area
8              Maryland       4                                 Horn Point
9              Maryland       6                    Howard County Near Road
10             Maryland       9                                 Hagerstown
    annl_mn                   geometry
1  9.661741 POINT (-76.95807 38.89557)
2  8.073839 POINT (-77.01318 38.92185)
3  8.695211 POINT (-76.95343 38.89477)
4  8.934325 POINT (-77.01282 38.87516)
5  6.922022 POINT (-75.79732 39.30502)
6  6.403240   POINT (-79.012 39.70595)
7  7.459565 POINT (-75.86477 39.70298)
8  6.918579 POINT (-76.14101 38.58752)
9  8.705889 POINT (-76.84611 39.14313)
10 7.068597 POINT (-77.72024 39.56418)

DMV Air Quality

ggplot() +
  geom_sf(data = dmv_county, color = "grey25", fill = "grey75", alpha = 0.55) +
  geom_sf(data = dmv_air, aes(colour = annl_mn), size = 2) + 
  scale_color_gradient(low = "green", high = "red") +
  labs(color = "Average PM2.5", x = "", y = "") +
  theme_bw()

Spatial Interpolation: Closest Observation

nw <- data.frame(long = -77.0654, lat = 38.9097)
nw <- st_as_sf(nw, coords = c("long", "lat"), crs = crs(dmv_county))
(nw$pm25_nav <- dmv_air$annl_mn[which.min(st_distance(nw, dmv_air))])

[1] 8.073839

Spatial Interpolation: IDW

(nw$pm25_idw <- sum(dmv_air$annl_mn / st_distance(nw, dmv_air)^2) / 
                          sum(1 / st_distance(nw, dmv_air)^2))

8.845472 [1]

library(gstat)
gstat_idw <- gstat(id = "pm25_idw", formula = annl_mn ~ 1, locations = dmv_air)
predict(gstat_idw, newdata = nw)

[inverse distance weighted interpolation]

Simple feature collection with 1 feature and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -77.0654 ymin: 38.9097 xmax: -77.0654 ymax: 38.9097
Geodetic CRS:  NAD83
  pm25_idw.pred pm25_idw.var                 geometry
1      8.845469           NA POINT (-77.0654 38.9097)

Spatial Interpolation: KNN

gstat_knn <- gstat(id = "pm25_knn", formula = annl_mn ~ 1, locations = dmv_air, nmax = 3,
             set = list(idp = 0)) # idp: inverse distance power; idp=0 means all three neighbors are equally weighted.
(nw$pm25_knn <- predict(gstat_knn, newdata = nw)$pm25_knn.pred)

[inverse distance weighted interpolation]

[1] 8.886756

Spatial Interpolation

Create locations for prediction

dmv_bbox <- st_bbox(dmv_county)
x <- seq(dmv_bbox$xmin, dmv_bbox$xmax, length.out = 50)
y <- seq(dmv_bbox$ymin, dmv_bbox$ymax, length.out = 25)
dc_locs <- expand.grid(x = x, y = y)
dc_pred0 <- st_as_sf(dc_locs, coords = c("x", "y"), crs = crs(dmv_county))
dc_pred <- st_intersection(dmv_county, dc_pred0)

Spatial Interpolation

ggplot() +
  geom_sf(data = dmv_county, color = "grey25", fill = "grey75", alpha = 0.55) +
  geom_sf(data = dc_pred, color = "blue", size = 1, shape = 4) + 
  labs(x = "", y = "") +
  theme_bw()

Spatial Interpolation

Using IDW

dc_pred$pm25_idw <- predict(gstat_idw, newdata = dc_pred)$pm25_idw.pred

[inverse distance weighted interpolation]

Using KNN

dc_pred$pm25_knn <- predict(gstat_knn, newdata = dc_pred)$pm25_knn.pred

[inverse distance weighted interpolation]

Spatial Interpolation

ggplot() +
  geom_sf(data = dmv_county, color = "grey25", fill = "grey75", alpha = 0.55) +
  geom_sf(data = dc_pred, aes(color = pm25_idw), size = 1.5) + 
  geom_sf(data = dmv_air, aes(color = annl_mn), size = 3) +
  scale_color_gradient(low = "green", high = "red") +
  labs(color = "PM2.5", x = "", y = "") +
  theme_bw()

Spatial Interpolation with Raster Data

library(raster)
dmv_raster <- raster(ext = extent(dmv_bbox), res = c(0.1, 0.1))
dmv_raster_sf <- st_as_sf(rasterToPolygons(dmv_raster)) %>%
  st_transform(crs = crs(dmv_county)) %>%
  st_intersection(dmv_county)

Raster Data

ggplot() +
  geom_sf(data = dmv_county, color = "grey35", fill = "grey75", alpha = 0.55) +
  geom_sf(data = dmv_raster_sf, color = "blue", size = 0.5) + 
  labs(x = "", y = "") +
  theme_bw()

Spatial Interpolation with Raster Data

dmv_raster_sf$pm25_idw <- predict(gstat_idw, newdata = dmv_raster_sf)$pm25_idw.pred

[inverse distance weighted interpolation]

Code

ggplot() +
      geom_sf(data = dmv_county, color = "grey25", fill = "grey75", alpha = 0.55) +
      geom_sf(data = dmv_raster_sf, aes(fill = pm25_idw), size = 1) +
      geom_sf_label(data = dmv_air, aes(label = round(annl_mn, 1), fill = annl_mn), size = 5, alpha=0.75) +
      scale_color_gradient(low = "green", high = "red") +
      scale_fill_gradient(low = "green", high = "red") +
      labs(fill = "PM2.5", x = "", y = "") +
      guides(color = guide_legend(order = 1)) +
      theme_bw() +
      theme(legend.title = element_text(size = 12),
            legend.text = element_text(size = 10))

Kriging

Kriging is an interpolation method named after the South African mining engineer, Danie G. Krige.
It’s originated in the areas of mining and geo-statistics that involve spatially correlated data.
Kriging assumes that the spatial variation of a variable is composed of three components: trend (structural component), spatially correlated variation, noise (random error)
Ordinary kriging: relies purely on a spatial error process to make predictions.
Universial kriging: allows spatial trend terms and even location-specific covariates to shape the prediction

Kriging是一种先进的空间插值方法，以南非采矿工程师Danie G. Krige的名字命名。Kriging最早用于采矿和地质统计领域，特别适合分析具有空间相关性的地理数据。与简单的距离或邻近方法不同，Kriging不仅考虑到观测点的位置，还将空间变化分解为多个组成部分，以生成更加精确的预测。

Kriging 的来源和应用 Kriging的概念起源于采矿和地质统计学，主要用于处理具有空间自相关的数据。在这些领域中，了解地理空间中变量的变化对于矿产勘探和资源评估至关重要。Kriging现已广泛应用于环境科学、土壤研究、气象学等多个学科中。
Kriging 的三大组成部分 Kriging假定变量的空间变化由三个主要部分组成：

趋势（Trend or Structural Component）：代表变量的整体趋势或结构性变化，例如气温随着纬度的变化趋势。空间相关变化（Spatially Correlated Variation）：捕捉局部的空间依赖性，即相邻位置之间的相似性。这是Kriging的核心部分，反映了变量的空间自相关。噪声（Noise or Random Error）：表示不可预测的随机误差或本地波动。噪声成分帮助我们识别预测中的不确定性。通过将空间变化分解为这三部分，Kriging能够生成更加准确和可靠的插值结果。

Kriging 的主要类型 Kriging 方法有多种变体，其中普通Kriging和泛Kriging是最常用的两种类型。

普通Kriging (Ordinary Kriging)：普通Kriging假设空间变异是随机的，仅依赖空间误差过程来生成预测结果。它不考虑特定位置的趋势，只关注位置之间的空间自相关。

泛Kriging (Universal Kriging)：泛Kriging不仅依赖于空间误差过程，还允许引入空间趋势项和特定位置的协变量。这意味着泛Kriging可以适应复杂的空间模式，比如在土壤研究中同时考虑地形和植被等影响因素。

这两种方法各有其适用场景：普通Kriging适合用于空间结构简单的情况，而泛Kriging则适合更复杂的地理模式。

Ordinary Kriging

Suppose we want to predict y(s_0) based observed data y(s_1), y(s_2), ..., y(s_n). The ordinary kriging estimator is defined as the linear unbiased estimator: \hat{y}(s_0) = \sum_{i=1}^n \lambda_i y(\mathbf{s}_i)

that minizes the mean squared prediction error defined as E[(\hat{y}(s_0) - y(s_0))^2]

Ordinary Kriging 是一种空间插值方法，用于预测未观测位置的变量值。它基于已有观测数据，通过加权平均估算目标位置的值，并且在计算中确保预测的无偏性和最小均方误差。接下来，我们深入了解普通Kriging的预测公式及其原理。

普通Kriging的预测公式假设我们想根据观测数据y(s_1), y(s_2), ..., y(s_n)来预测位置s_0的值。普通Kriging的预测公式如下：

在这个公式中，\hat{y}(s_0) 表示对位置s_0的预测值。 y(\mathbf{s}_i) 是位置s_i的观测值。 \lambda_i 是每个观测值的加权系数，根据位置s_i和目标位置s_0之间的空间相关性来确定。这些权重\lambda_i是关键，它们使得Kriging预测能够更好地反映空间结构。通过调整这些权重，我们可以让预测值对相邻观测点的值更敏感，捕捉空间自相关性。

最小化均方预测误差普通Kriging不仅关注预测准确性，还通过最小化均方预测误差来优化预测值。这一误差定义如下：

目标是使\hat{y}(s_0)的均方误差最小，从而在统计上获得最优的预测。通过最小化预测误差，普通Kriging确保了预测值与真实值之间的差异尽可能小，这在空间数据不完全或数据噪声较大时尤为重要。在实际计算中，普通Kriging通过方差和协方差矩阵来构建\lambda_i的最优组合，使得预测无偏且具有最小的均方误差。

Universial Kriging

y(\mathbf{s}) = \mu(\mathbf{s}) + \omega(\mathbf{s}) + \epsilon(\mathbf{s})

\mu(\mathbf{s}) = \mathbf{X}(\mathbf{s})\beta

\omega(\mathbf{s}) \sim N(0, \sigma^2\mathbf{H}(\theta))

Partial sill: \sigma^2 is the variance driven by distance between two observations
Decay: \theta is the level of correlation among observations based on distance.

\sigma(\mathbf{s}) \sim N(0, \tau^2\mathbf{I})

Nugget: \tau^2 is the error variance that is independent from spatial separation

Universal Kriging 是一种更为灵活的Kriging方法，适用于空间数据中包含明显趋势的情况。与普通Kriging不同，Universal Kriging不仅考虑空间自相关，还允许引入空间趋势和位置特定的协变量，以更精确地捕捉空间变化。Universal Kriging的预测模型可以分解为三个主要组成部分：趋势、空间相关性和随机误差。

Universal Kriging 的模型结构 Universal Kriging 的基本模型可以表示为：y(s)=μ(s)+ω(s)+ϵ(s)

\mu(\mathbf{s})：表示位置s的趋势成分。 \omega(\mathbf{s})：表示空间相关性成分。 \epsilon(\mathbf{s})：表示独立随机误差成分（即噪声）。这种模型结构使Universal Kriging能够适应更复杂的空间模式，特别适合存在空间趋势的情况，例如地形影响气候的变化或植被影响土壤湿度的情况。

趋势成分 \mu(\mathbf{s}) 在Universal Kriging中，趋势成分\mu(\mathbf{s})通常被定义为观测点的协变量\mathbf{X}(\mathbf{s})与系数\beta的线性组合：

\mathbf{X}(\mathbf{s}) 可以是位置特定的特征或协变量（例如土壤类型、坡度、气温等）。 \beta 是协变量的回归系数，反映每个特征对预测结果的影响。通过引入协变量，Universal Kriging能够捕捉空间上的显著趋势，为预测提供更多信息。

空间相关性成分 \omega(\mathbf{s}) 空间相关性成分\omega(\mathbf{s})建模为零均值的正态分布，协方差矩阵由空间距离驱动：

部分基台 (Partial Sill)：\sigma^2表示空间变化的方差，反映了相邻观测点间随距离增加而逐渐减弱的相关性。衰减 (Decay)：\theta表示距离驱动的相关性，控制相邻位置之间的相关性衰减速度。较小的\theta值表示更快的衰减，较大的值表示更缓慢的空间衰减。这个成分捕捉到变量的空间自相关性，使得Universal Kriging能够处理位置间的空间依赖。

噪声成分 \epsilon(\mathbf{s}) 噪声成分\epsilon(\mathbf{s})表示随机误差，假定为零均值正态分布：

基块 (Nugget)：\tau^2表示独立于空间位置的误差方差。Nugget成分用于捕捉小尺度波动或测量误差，即即使在零距离的观测点也可能存在不完全相等的值。 Nugget项反映了不可解释的随机误差，有助于处理小尺度随机波动。

Variogram

In order to make spatial predictions, kriging goes through two steps:
- Training: using the variograms and covariance functions to estimate the spatial autocorrelation
- Prediction: predicting the unknown values

在Kriging插值方法中，半变异函数（Variogram）是一个关键工具，用于量化空间自相关性。Kriging的预测过程可以分为两个步骤：训练和预测。在训练阶段，半变异函数帮助我们估计数据中的空间自相关程度，从而构建适合的空间预测模型。接下来，我们将详细介绍这两个步骤。

Kriging的两个步骤 Kriging插值过程主要包括以下两个步骤：

训练 (Training)：通过半变异函数（Variogram）和协方差函数（Covariance Function）来估计数据的空间自相关性。这一步帮助我们理解观测点之间的关系，建立一个能够捕捉空间依赖性的模型。

预测 (Prediction)：使用训练阶段建立的模型，对未知位置进行预测。基于空间自相关性，我们可以根据邻近位置的数据估算未观测位置的值。

半变异函数 (Variogram) 的作用半变异函数是Kriging中用于量化空间自相关性的关键工具。它描述了随着空间距离增加，变量值的相似性如何变化。具体来说，半变异函数反映了两个位置之间的距离与变量差异之间的关系。

原理：当两个位置的距离较近时，半变异值较小，表示它们的属性值更相似；随着距离增加，半变异值增大，表示它们的属性值差异增加。结构：半变异函数通常包含三个关键参数：基块（Nugget）、部分基台（Partial Sill）和范围（Range）。这些参数帮助我们描述变量的空间自相关模式。通过半变异函数，我们可以了解数据在不同空间距离下的相似程度，从而选择适合的Kriging模型。

训练与预测的关系在训练阶段，半变异函数提供了关于空间结构的统计信息，使我们能够为每个观测点分配适当的权重。在预测阶段，Kriging利用这些权重来进行加权平均，从而在未知位置生成合理的预测值。

这种训练和预测的过程使Kriging能够根据已有数据的空间依赖性来进行准确的插值预测。

Variogram

Theoretical variogram:

Variogram

To obtain the theoretical variogram, we can use empirical semivariogram

Code

library(gstat)
boston_evgm <- variogram(MEDV ~ 1, boston)
plot(boston_evgm, type = 'o', main = "Variogram of Boston Median House Values")

Variogram

Code

v_init <- vgm(psill = 50, model = "Sph",
              range = 5, nugget = 50)
plot(boston_evgm, v_init, cex = 1.5)

Variogram

Then, we can choose a theoretical variogram model that best fits the empirical semivariogram.
- Gaussian, Exponential, and Spherical models
- Model parameters (nugget, sill, range) are estimated by fitting the chosen model to the empirical semivariogram.

Variogram

fittted_vgm <- fit.variogram(object = boston_evgm,
                    model = vgm(psill = 50, model = "Sph",
                                range = 5, nugget = 50))
fittted_vgm

  model    psill    range
1   Nug  0.00000 0.000000
2   Sph 93.13453 3.365808

Variogram

plot(boston_evgm, fittted_vgm, cex = 1.5)

Kriging

kmod <- gstat(formula = MEDV ~ 1, data = boston, model = fittted_vgm)

pred_loc <- st_as_sf(data.frame(long = -71.1456, lat = 42.3915), 
                       coords = c("long", "lat"), crs = crs(boston))

predict(kmod, pred_loc)

[using ordinary kriging]

Simple feature collection with 1 feature and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -71.1456 ymin: 42.3915 xmax: -71.1456 ymax: 42.3915
Geodetic CRS:  NAD27
  var1.pred var1.var                 geometry
1  21.42989 18.91051 POINT (-71.1456 42.3915)

Spatial Regression Model

OLS
- Gauss-Markov assumption: Cov(\epsilon_i, \epsilon_j) = 0, \forall i \neq j
Spatial autocorrelation
- Spatially autocorrelated errors violates the OLS assumption, indicating that the model is misspecified.

在地理空间数据中，普通最小二乘法（Ordinary Least Squares, OLS）是基础的回归方法。然而，空间数据的特点会引入空间自相关，使得OLS模型不再满足经典的Gauss-Markov假设。这种空间自相关性意味着标准OLS模型可能无法捕捉到变量之间的空间关系，从而导致模型设定错误。接下来，我们详细探讨这一问题。

OLS与Gauss-Markov假设在传统的OLS回归模型中，Gauss-Markov假设要求误差项之间相互独立，且满足以下条件：含义：这个假设意味着每个观测点的误差项\epsilon_i都不应该与其他观测点的误差项\epsilon_j相关。换句话说，观测点之间的误差没有空间依赖性。结果：在这个假设成立的前提下，OLS回归是最佳线性无偏估计（BLUE），能够生成精确且无偏的估计结果。
空间自相关性与OLS假设的冲突然而，在地理空间数据中，观测点之间往往存在空间自相关性。当数据呈现空间聚集或空间依赖时，误差项\epsilon_i和\epsilon_j之间不再独立，这违背了Gauss-Markov假设。

空间自相关的影响：如果模型中未考虑空间自相关性，OLS估计结果可能会产生偏差，导致模型误设。空间自相关性表明观测点之间存在某种关系，而OLS假设忽略了这种空间依赖性。模型设定错误：当空间自相关性未被考虑时，OLS模型的残差可能表现出空间模式。这种情况通常表示模型遗漏了重要的空间因素，导致回归结果不准确。在这种情况下，我们需要采用空间回归模型，如空间滞后模型（Spatial Lag Model）或空间误差模型（Spatial Error Model），以充分考虑地理数据中的空间依赖性，从而提高模型的解释力和预测力。

OLS and Spatial Errors

lm_fit <- lm(log(MEDV) ~ RM + CRIM, data = boston)
summary(lm_fit)


Call:
lm(formula = log(MEDV) ~ RM + CRIM, data = boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.03618 -0.11069  0.01723  0.13889  1.31727 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.121864   0.112901   9.937   <2e-16 ***
RM           0.315514   0.017661  17.865   <2e-16 ***
CRIM        -0.019438   0.001443 -13.474   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2721 on 503 degrees of freedom
Multiple R-squared:  0.5587,    Adjusted R-squared:  0.557 
F-statistic: 318.4 on 2 and 503 DF,  p-value: < 2.2e-16

Spatially Autocorrelated Errors

boston$residuals <- residuals(lm_fit)

ggplot(boston, aes(fill = residuals)) +
  geom_sf() +
  scale_fill_steps(low = "grey", high = "brown") + 
  theme_bw()

Spatially Autocorrelated Errors

nb <- poly2nb(boston) 
nbw <- nb2listw(nb, style = "W")
moran_res <- moran.mc(boston$residuals, nbw, 1000)
moran_res


    Monte-Carlo simulation of Moran I

data:  boston$residuals 
weights: nbw  
number of simulations + 1: 1001 

statistic = 0.54595, observed rank = 1001, p-value = 0.000999
alternative hypothesis: greater

Spatial Lag Model

Also known as the spatial autoregressive model (SAR)
Incorporating spatial dependence in the outcome y = \rho \mathbf{W} y + \mathbf{X} \beta + \epsilon
- \rho: Spatial autoregressive coefficient
- \mathbf{W}: Spatial weights matrix

在空间回归分析中，空间滞后模型（Spatial Lag Model），也称为空间自回归模型（Spatial Autoregressive Model, SAR），是用于处理因变量中空间依赖性的一种模型。空间滞后模型通过在回归方程中引入因变量的空间滞后项，捕捉到空间单元之间的相互影响。

空间滞后模型的方程空间滞后模型的基本形式如下：

y: 因变量向量，代表每个空间单元的观测值。 \rho: 空间自回归系数，衡量因变量的空间滞后项对因变量的影响强度。 \mathbf{W}: 空间权重矩阵，用于定义空间单元之间的相邻关系和空间依赖性。 \mathbf{X}: 自变量矩阵，包括每个空间单元的特征。 \beta: 回归系数，表示自变量对因变量的影响。 \epsilon: 误差项，假定独立且无空间相关性。 2. 解释空间滞后项 \rho \mathbf{W} y 空间滞后项 \rho \mathbf{W} y 表示每个空间单元的观测值如何受到其邻近单元观测值的影响。

空间自回归系数 \rho：如果\rho为正，表示相邻单元之间的观测值具有正空间自相关性，即高值倾向于与高值相邻，低值倾向于与低值相邻；如果\rho为负，则表示负空间自相关性。空间权重矩阵 \mathbf{W}：空间权重矩阵定义了每个单元与其邻居的关系。不同的权重矩阵（如邻接矩阵或反距离矩阵）可以反映不同的空间依赖结构。通过引入空间滞后项，空间滞后模型能够捕捉到因变量中的空间依赖性，使模型更好地反映空间单元之间的相互影响。

空间滞后模型的应用空间滞后模型适用于因变量表现出显著空间聚集效应的场景：

房地产市场：房价往往受到相邻区域房价的影响，因此在房价预测中可以使用空间滞后模型。社会经济研究：人均收入、失业率等变量通常表现出空间聚集性，使用空间滞后模型可以更好地分析这些变量的空间分布。公共健康：疾病传播可能受到相邻地区病例数的影响，空间滞后模型可以帮助分析疾病的空间扩散模式。

Spatial Lag Model

library(spatialreg)
lagsar_fit <- lagsarlm(log(MEDV) ~ RM + CRIM, data=boston, nbw)
summary(lagsar_fit)


Call:lagsarlm(formula = log(MEDV) ~ RM + CRIM, data = boston, listw = nbw)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.6189521 -0.1004647 -0.0049353  0.0898597  1.0161974 

Type: lag 
Coefficients: (asymptotic standard errors) 
              Estimate Std. Error z value  Pr(>|z|)
(Intercept) -0.2810648  0.0838293 -3.3528 0.0007999
RM           0.1977215  0.0132364 14.9377 < 2.2e-16
CRIM        -0.0070405  0.0010320 -6.8219 8.986e-12

Rho: 0.68944, LR test value: 384.41, p-value: < 2.22e-16
Asymptotic standard error: 0.027224
    z-value: 25.325, p-value: < 2.22e-16
Wald statistic: 641.35, p-value: < 2.22e-16

Log likelihood: 134.3822 for lag model
ML residual variance (sigma squared): 0.030777, (sigma: 0.17543)
Number of observations: 506 
Number of parameters estimated: 5 
AIC: -258.76, (AIC for lm: 123.65)
LM test for residual autocorrelation
test value: 6.7847, p-value: 0.0091942

Spatial Lag Model

boston$residuals2 <- residuals(lagsar_fit)
moran.mc(boston$residuals2, nbw, 1000)


    Monte-Carlo simulation of Moran I

data:  boston$residuals2 
weights: nbw  
number of simulations + 1: 1001 

statistic = 0.054806, observed rank = 988, p-value = 0.01299
alternative hypothesis: greater

Spatial Error Model

Also called spatial moving average models
Incorporating spatial dependence in the error term y = \mathbf{X} \beta + u u = \lambda \mathbf{W} u + \epsilon
- u: Spatially correlated error term
- \lambda: Spatial autoregressive coefficient for the error term
- \mathbf{W}: Spatial weights matrix

在空间回归分析中，空间误差模型（Spatial Error Model），也称为空间移动平均模型（Spatial Moving Average Model），用于处理模型中误差项的空间依赖性。与空间滞后模型不同，空间误差模型并不在因变量中引入空间依赖，而是通过误差项的结构来捕捉空间相关性。

空间误差模型的方程空间误差模型的基本形式可以分为两部分：

回归方程：空间依赖的误差项： λWu+ϵ

y: 因变量向量，表示每个空间单元的观测值。 \mathbf{X}: 自变量矩阵，包括每个空间单元的特征。 \beta: 回归系数，表示自变量对因变量的影响。 u: 空间相关的误差项，反映了未解释的空间依赖性。 \lambda: 误差项的空间自回归系数，控制空间误差项的强度。 \mathbf{W}: 空间权重矩阵，用于定义空间单元之间的相邻关系。 \epsilon: 随机误差项，假定独立且无空间自相关性。 2. 解释误差项 u = \lambda \mathbf{W} u + \epsilon 在空间误差模型中，误差项u本身存在空间相关性，由空间权重矩阵\mathbf{W}定义。

空间自回归系数 \lambda：该系数表示误差项中空间相关性的强度。如果\lambda为正，表示误差项具有正空间自相关性，即邻近单元的误差项之间的值更倾向于相似；如果\lambda为负，表示误差项存在负空间自相关性。空间权重矩阵 \mathbf{W}：权重矩阵定义了每个单元与相邻单元的关系，这一结构使得误差项可以反映邻近位置的未解释空间依赖。通过这种方式，空间误差模型捕捉到了模型中未被解释的空间依赖，从而提高了估计的准确性。

Spatial Error Model

ersar_fit <- errorsarlm(log(MEDV) ~ RM + CRIM, data=boston, nbw, tol.solve=1.0e-30)
summary(ersar_fit)


Call:errorsarlm(formula = log(MEDV) ~ RM + CRIM, data = boston, listw = nbw, 
    tol.solve = 1e-30)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.9022086 -0.0739408 -0.0003695  0.0841594  0.9155054 

Type: error 
Coefficients: (asymptotic standard errors) 
              Estimate Std. Error z value  Pr(>|z|)
(Intercept)  1.4309386  0.1011931 14.1407 < 2.2e-16
RM           0.2544088  0.0141162 18.0224 < 2.2e-16
CRIM        -0.0071259  0.0012310 -5.7886 7.099e-09

Lambda: 0.85125, LR test value: 390.82, p-value: < 2.22e-16
Asymptotic standard error: 0.02544
    z-value: 33.462, p-value: < 2.22e-16
Wald statistic: 1119.7, p-value: < 2.22e-16

Log likelihood: 137.5884 for error model
ML residual variance (sigma squared): 0.027826, (sigma: 0.16681)
Number of observations: 506 
Number of parameters estimated: 5 
AIC: -265.18, (AIC for lm: 123.65)

Spatial Error Model

boston$residuals3 <- residuals(ersar_fit)
moran.mc(boston$residuals3, nbw, 1000)


    Monte-Carlo simulation of Moran I

data:  boston$residuals3 
weights: nbw  
number of simulations + 1: 1001 

statistic = -0.036993, observed rank = 93, p-value = 0.9071
alternative hypothesis: greater

Spatial Error Model

ggplot(boston, aes(fill = residuals3)) +
  geom_sf() +
  scale_fill_steps(low = "grey", high = "brown") + 
  theme_bw()

Geographically Weighted Regression (GWR)

Allowing the regression coefficients to vary spatially y_i = \beta_0(u_i, v_i) + \sum_{k=1}^p \beta_k(u_i, v_i) x_{ik} + \epsilon_i
- (u_i, v_i): Spatial coordinates of location i
- \beta_k(u_i, v_i): Spatially varying regression coefficients
- x_{ik}: Value of the k-th explanatory variable at location i
- \epsilon_i: Error term at location i
Estimating local regression coefficients at each location
Capturing spatial non-stationarity and local variations in relationships

Geographically Weighted Regression (GWR) 是一种空间回归模型，它允许回归系数在空间上变化，从而捕捉不同位置之间的局部关系。与传统的回归模型不同，GWR假设自变量对因变量的影响可能因地理位置而异，因此可以有效捕捉空间数据中的空间非平稳性（Spatial Non-Stationarity）。

GWR 模型的方程 GWR的基本方程如下： (u_i, v_i)：位置 i 的空间坐标，用于识别地理位置。 \beta_k(u_i, v_i)：在位置 i 处第 k 个解释变量的空间回归系数。GWR允许这些系数随空间位置而变化，体现出不同区域的关系模式。 x_{ik}：在位置 i 处第 k 个解释变量的观测值。 \epsilon_i：位置 i 的误差项。 GWR 通过允许回归系数在空间上变化，能够捕捉到因变量与自变量之间的局部关系，而非假设这种关系在整个空间上是固定的。
GWR 的特点与作用 GWR的核心在于估计每个位置的局部回归系数，使得模型能够适应不同位置的空间特征：

捕捉空间非平稳性：传统回归模型假设系数在空间上是固定的，而GWR放宽了这一假设，允许不同位置的系数变化，从而更好地反映不同区域的局部关系。本地化变化：GWR能够揭示自变量对因变量的影响如何在不同区域内发生变化，这在分析具有高度空间异质性的现象时尤其有用。例如，在房价分析中，土地面积、人口密度、教育水平等变量对房价的影响可能会随地区而变化。GWR可以估计每个区域的回归系数，从而捕捉这些差异。

What We Didn’t Get To Cover

Point process
Raster data
Hierarchical models for spatial data
Machine learning for spatial data

点过程 (Point Process) 点过程用于分析空间中离散事件的分布，尤其适合处理事件在特定区域中的发生位置。

应用：点过程广泛应用于生态学（动物栖息地分布）、犯罪分析（犯罪事件分布）以及流行病学（疾病爆发的空间分布）等领域。目标：点过程模型可以帮助我们理解事件发生的空间模式，如聚集、均匀或随机分布，进而揭示潜在的影响因素。 2. 栅格数据 (Raster Data) 栅格数据特别适合处理连续的空间变量，通常用于表示地理表面的特征，如高程、温度和土地覆盖。

特点：栅格数据将空间划分为一系列网格单元，每个单元包含变量的数值。这种数据结构便于处理和分析大规模、连续变化的地理数据。应用：栅格数据在遥感分析、环境科学、土地利用等领域中应用广泛，适合生成空间分布图和进行地表分析。 3. (Bayesian) Hierarchical Models for Spatial Data 层次贝叶斯模型提供了一种灵活的方法来处理空间数据中的多层次和多尺度问题。通过贝叶斯框架，我们可以将空间数据建模为多个层次的结构，捕捉到不同层级的空间变化。

推荐书籍：Banerjee, S., Carlin, B.P., & Gelfand, A.E. (2014). Hierarchical Modeling and Analysis for Spatial Data (2nd ed.). 该书详细介绍了空间数据的层次建模方法，特别适合想深入学习贝叶斯空间模型的读者。应用：层次贝叶斯模型广泛用于生态学、环境监测和公共卫生等领域，能够处理空间数据中复杂的依赖结构和不确定性。 4. 空间数据的机器学习 (Machine Learning for Spatial Data) 随着数据科学的发展，机器学习方法在空间数据分析中发挥了越来越重要的作用。机器学习能够处理大规模、复杂的空间数据，并且适用于非线性和高维度的分析任务。

应用：机器学习在土地分类、图像识别、预测分析等方面表现突出。例如，深度学习可以用于卫星图像分析，而随机森林可以处理大规模的地理空间数据。优势：机器学习方法适合自动化处理和模式识别，为空间数据分析带来了新的可能性和应用场景。