1. Overview of spatial data problems

본 내용은 아래의 책을 참조하여 작성하였습니다

https://www.researchgate.net/publication/224839687_Hierarchical_Modeling_and_Analysis_of_Spatial_Data

---

1.1 Introduction to spatial data and models

 

• Spatial data는 크게 3가지 type으로 구분된다.

1. point-referenced data, where $Y(s)$ is a random vector at a location $s \in \mathcal{R}^r$, where $s$ varies continuously over $D$, a fixed subset of $\mathcal{R}^r$ that contains an $r$-dimensional rectangle of positive volume
2. areal data, where $D$ is again a fixed subset (of regular or irregular shape), but now partitioned into a finite number of areal units with well-defined boundaries
3. point pattern data, where now $D$ is itself random; its index set gives the locations of random events that are the spatial point pattern. $Y(s)$ itself can simply equal 1 for all $s \in D$ (indicating occurrence of the event), or possibly give some additional covariate information (producing a marked point pattern process).

 

• 필자는 간단하게 다음과 같이 이해했다.

1. point-referenced data : 위치를 좌표로 나타낸 것 (아마도 가장 많이 만나게 될 듯)
2. areal data : 위치가 특정 기준으로 나눠져 있는 것
3. point pattern data : 위치가 random인 것

 

• 또한 각 type을 뜻하는 용어가 여러가지 있다.

1. point-referenced data, point-level data, geocoded, geostatistical
2. areal data, lattice
3. (spatial) point pattern data, (spatial) point process

---

1.1.1 Point-level models

 

• point-level data를 다루기 위해서는 다음과 같은 가정을 필요로 한다.

The covariance between the random variables at two locations depends on the distance between the locations.

 

 

• 이런 covariance를 위해 주로 사용되는 것이 exponential model이다. 즉, $$Cov(Y(s_i), Y(s_j)) \equiv C(d_{ij}) = \sigma^2 e^{-\phi d_{ij}} \text{ for } i \neq j.$$ 여기서 $d_{ij}$는 $s_i$와 $s_j$ 간의 거리, $\sigma^2$와 $\phi$는 양의 값을 가지는 partial sill과 decay parameter이다.  
• $C(d_{ij}) = Var(Y(s_i))$는 흔히 $\tau^2 + \sigma^2$로 확장된다. 여기서 $\tau^2 > 0$는 nugget effect, $\tau^2 + \sigma^2$는 sill이라 부른다.

 

• covariance와 distance를 나타낸 plot은 covariogram이라 부른다.

 

• 이런 variancce와 covariance의  assumptions에 joint distributional model을 적용하면 일반적인 likelihood inference가 가능해진다.
• 가장 편한 approach는 multivariate normal을 가정하는 것이다. 즉, $\mathbf{Y} \equiv {Y(s_i)}$ at known locations $s_i, i = 1, ..., n$일 때, $$\mathbf{Y} | \mu, \Theta \sim N_n(\mu 1, \Sigma(\Theta)).$$ 여기서 $N_n$은 $n$-dimensional normal distribution, $\mu$는 (constant) mean level, $1$은 vector of ones, $(\Sigma(\Theta))_{ij}$는 $Y(s_i)$와 $Y(s_j)$간의 covariance, $\Theta = (\tau^2, \sigma^2, \phi)^T$이다.
• $\Sigma$의 형태를 고려함에 있어 가장 쉬운 선택은 isotropic covariance functions를 따르는 것이다. isotropic covariance functions는 spatial correlation이 오로지 $s_i$와 $s_j$ 간의 거리에 의존한다는 가정을 한 함수이다. 즉, $$(\Sigma(\theta))_{ij} = \sigma^2 \text{exp}(-\phi d_{ij}) + \tau^2 I(i = j), \sigma^2 > 0, \phi > 0, \tau^2 >0$$ 와 같은 형태이다.
• 다른 선택으로 위 식에서 $d$에 power를 준 형태, the spherical, the Gaussian, the Matern 등이 있다.

---

1.1.2 Areal models

 

• Areal data에서는 zip codes, counties, etc를 나타내기 위한 geographic regions (or blocks)를 $B_i$로 나타낸다. 또한 data는 block 내에서 variable의 합이나 평균이 된다.

 

• spatial association을 설명하기 위해 map 상의  blocks의 배열을 기반으로 한 neighborhood structure를 정의해야한다. 이런 neighborhood structure를 정의했다면 model은 autoregressive time series model과 비슷하게 고려할 수 있다.
• neighborhood information을 포함한 두 가지 가장 유명한 model은 simultaneously and conditionally autoregressive model (SAR and CAR)이다.
• SAR model은 likelihood methods를 사용하기에 계산적으로 편리함이 있다.
• 반면에 CAR model은 Bayesian model fitting과 결합한 Gibbs sampling을 사용하기에 계산적으로 편리함이 있다. 또한 주로 vector of spatially varying random effects $\Phi = (\phi_1, ..., \phi_n)^T$를 사용한다.
• 예를 들어, $Y_i \equiv Y(B_i)$라 하고 $Y_i \sim N(\phi_i, \sigma^2)$일 때, CAR model에서는 $$\phi_i | \Phi_{(-i)} \sim N(\mu + \sum^n_{k = 1} a_{ik} (\phi_k - \mu), \tau^2_i)$$를 고려한다. 여기서 $\Phi_{(-i)} = \{ \phi_k : k \neq i \}$, $\tau^2_i$는 conditional variance, $a_{ik}$는 known (or unknown) constant such that $a_{ii} = 0$.
• $A = (a_{ik}), M = Diag(\tau^2_1, ..., \tau^2_n)$라 두자. Brook's Lemma에 의해, $$p(\Phi) \propto \text{exp} \{-(\Phi - \mu 1)^T M^{-1} (I - A) (\Phi - \mu 1) / 2 \}$$를 보일 수 있다. 여기서 $1$은 $n$-vector of 1's, $I$는 $n \times n$ identity matrix이다.
• 일반적으로, $A$와 $M$은 다음과 같이 계산한다. $$A = \rho ~ Diag(1/w_{i+})W \text{   and   } M^{-1} = \tau^{-2} ~Diag(w_{i+}).$$ 여기서 $\rho$는 spatial correlation parameter, $W = (w_{ij})$는 neighborhood matrix for the areal unit이다. 특히 $w_{ij}$는 다음과 같이 정의한다. $$ \begin{align} w_{ij} = \Bigg\{
          \begin{array}{ll}
         1, & \text{if  subregions } i \text{ and } j \text{ share a common boundary}, i \neq j \\
          0, & \text{otherwise}
          \end{array}
  \end{align} $$
• $Diag(w_{i+})$는 $w_{i+} = \sum_j w_{ij}$인 diagonal matrix이다.
• $\alpha \equiv (\rho, \tau^2)$라 두자. $\Phi$의 covariance matrix는 $$C(\alpha) = \tau^2 [Diag(w_{i+}) - \rho W]^{-1}$$ 로 계산된다.

---

1.1.3 Point process models

 

• Point process model에서 $Y(s)$는 일반적으로 1 (indicating occurrence of the pattern)이며 추가적인 covariate information을 사용하게 된다.

 

• 해당 model에서 이루어 지는 task는 특정 location에서 event의 occurrence에 대한 expected value, 또는 어떤 clustering의 결과에 대한 조사가 있다.

 

• 자주 사용되는 함수가 있는데 Ripley's $K$ function이라고 하며 $$K(d) = \frac{1}{\lambda} E[\text{number of points within } d \text{ of an arbitrary point}]$$ 이다. 여기서 $\lambda$는 intensity parameter (i.e. mean number of points per unit area)이다. $K$의 theoretical value는 특정 spatial point process model로 유명하다. 
• 모든 location에서 spatial dependence가 없다고 가정하면 $K(d) = \pi d^2$로 두어도 무방할 것이다. 단순히 반지름 $d$에 해당하는 원의 넓이에 arbitrary point의 $d$ 이내에 존재하는 point의 수가 비례하기 때문이며 이로 인해 식에서 $\lambda$가 나누어 지는 것이다.
• 일반적인 $K$의 estimator는 $$\hat K(d) = n^{-2} |A| \sum \sum^{i \neq j} p^{-1}_{ij} I_d(d_{ij})$$ 이다. 여기서 $n$은 영역 $A$ 내의 point의 수, $p_{ij}$는 중심이 $i$인 원의 $A$ 내의 $j$를 통과하는 비율(?), $I_d(d_{ij})$는 $d_{ij} < d$이면 1, 그렇지 않으면 0의 값을 갖는다.

---

1.2 Fundamentals of cartography

skip

---

1.3 Maps and geodesics in R

https://conservancy.umn.edu/handle/11299/200500

Supplemental Materials to Hierarchical Modeling and Analysis for Spatial Data, 2nd Edition

library(tidyverse)

# -------------------------------------------------------------------------

library(maps) # install.packages("maps")
detach("package:purrr") # for function map
mn.map <- map(database = "county", region = "minnesota")
mn.map %>% str

---

library(rgdal) # install.packages("maptools")
cities.shp <- readOGR(dsn = system.file("vectors", package = "rgdal"), layer = "cities")
plot(cities.shp)
cities.shp %>% str

---

library(maps)
library(mapdata) # install.packages("mapdata")
map(database = "world", regions = "Canada", xlim = c(-141, -53), ylim = c(40, 85), col = "gray90", fill = TRUE)

---

coloradoST <- read.table("HierarchicalModelingAndAnalysis_SupplementalFiles//ColoradoS-T.dat", header = TRUE)
library(RgoogleMaps) # install.packages("RgoogleMaps")
MyMap <- GetMap.bbox(lonR = range(coloradoST$Lon),
                     latR = range(coloradoST$Lat),
                     size=c(640,640), maptype = "hybrid")
PlotOnStaticMap(MyMap)

convert_points <- LatLon2XY.centered(MyMap, coloradoST$Lat, coloradoST$Lon)
points(convert_points$newX, convert_points$newY, col = "red", pch = 19)

---

library(sp)
library(rgdal)
SP_longlat <- SpatialPoints(coords = cbind(coloradoST$Lon, coloradoST$Lat), proj4string = CRS("+proj=longlat +datum=WGS84"))
SP_utm <- spTransform(SP_longlat, CRS("+proj=utm +zone=13 +datum=WGS84"))
plot(SP_utm)

---

library(fields) # install.packages("fields")
X1 <- coloradoST[1, c(1, 2)]
X2 <- coloradoST[2, c(1, 2)]
euclidean.dist = rdist(X1, X2)
spherical.dist = rdist.earth(X1,X2)

---

1.4 Exercises

skip