INLA는 잠재 가우스 모형 (latent gaussian model)에 특화된 방법이며 복잡한 사후 분포를 단순화된 라플라스 근사치를 이용하고 적분을 유한개의 가중 합으로 계산하기 때문에 MCMC에 비해 비교적 짧은 계산시간에 정확한 결과를 제공한다는 이점이 있다. 이러한 장점을 이용하여 2019년 서울에서 발생한 강간강제추행을 분석하기 위해 R-INLA (Integrated Nested Laplace Approximations)를 적용해보았다. 데이터는 서울 열린데이터광장에서 수집하였으며, BYM Model을 이용하여 공간 분석을 실시하였다.
BYM Model
서울 강간강제추행범죄는 count data이기 때문에 Poisson regression이 사용된다. \(i\)번째 자치구의 강간강제추행 발생 건수를 \(y_{i}\)라고 할 때, \[\begin{align*} y_{i} &\thicksim Poisson(\lambda_{i}), \;\; \lambda_{i} = E_{i}\rho_{i}\\ \eta_{i} &= \log{\rho_{i}} = b_{0} + u_{i} + v_{i} \end{align*}\]
- \(E_{i}\) : \(i\)번째 자치구에서 예상된 수
- \(\rho_{i}\) : \(i\)번째 자치구의 상대적 위험 (Relative risk)
- \(b_{0}\) : 모든 자치구에서의 평균 \(y\)
- \(u_{i}\) : 구조화된 공간 효과 (Structured spatial effect)
- \(v_{i}\) : 비구조화된 공간 효과 (Unstructured spatial effect)
Spatial Effect
BYM Model은 공간적으로 구조화된 요소(\(u_{i}\))와 비구조화된 요소(\(v_{i}\)) 두 가지를 모두 고려한다.
구조화된 공간 효과(\(u_{i}\))
\[\begin{align*} &u_{i} \vert u_{-i} \thicksim N(m_{i}, s^2_{i}),\\ &m_{i} = \frac{\sum_{j \in \mathcal{N}_{i}} u_{j}}{{\# \mathcal{N}_{i}}} ,\;\; s^2_{i} = \frac{\sigma^2_{u}}{\# \mathcal{N}_{i}}, \end{align*}\]
\(\mathcal{N}_{i}\) : \(i\) 번째 지역과 근접한 이웃 지역 집합
\(\# \mathcal{N}_{i}\) : \(i\) 번째 지역과 근접한 이웃 지역 갯수
- 근접한 이웃은 경계선이 맞닿아 있는 지역을 의미
\(m_{i}\) : \(i\) 번째 지역과 근접한 이웃 지역들의 \(u_{j}\) 평균
\(s^2_{i}\) : 근접한 이웃 수에 의존
- 이웃 수가 많으면 분산은 작아짐
- 강한 공간 상관관계가 존재하는 경우, 이웃을 많이 가지는 지역일수록 더 많은 정보를 가진다는 것
- 이웃 수가 많으면 분산은 작아짐
공간 상관관계는 iCAR을 사용
비구조화된 공간 효과(\(v_{i}\))
\[\begin{align*} v_{i} \thicksim N(0, \sigma^2_{v}) \end{align*}\]
- 순수한 과산포 (Overdispersion)에 대한 고려
BYM Model에서의 모수 (Parameter)는 \(log{\tau_{u}}\), \(log{\tau_{v}}\) , where \(\tau_{u} = 1/{\sigma^2_{u}}\), \(\tau_{v} = 1/{\sigma^2_{v}}\)
\(u_{i}\)와 \(v_{i}\)는 독립적으로 해석할 수 없으며, 오직 \(\xi_{i}=u_{i}+v_{i}\)로 공간 효과를 식별할 수 있음
Real Data
2019년 서울에서 발생한 강간강제추행범죄를 분석하기 위해 R-INLA (Integrated Nested Laplace Approximations)를 이용하여 BYM Model을 적용해보았다. 데이터는 서울 열린데이터광장에서 수집하였다.
Loading Data
pacman::p_load("maptools", # For readShapePoly
"spdep", # For poly2nb
"dplyr",
"RColorBrewer", # For brewer.pal
"INLA")
dat.2019 <- read.csv("2019_crime.csv",header=T)
# Convert rows in the order of ESPI_PK
dat.2019 <- dat.2019[order(dat.2019$ESRI_PK),]
head(dat.2019)
ESRI_PK year district rape pop_total pop_femal sec_fac
10 0 2019 도봉구 90 335631 171670 10
12 1 2019 은평구 194 484546 251186 22
6 2 2019 동대문구 168 363023 184533 21
20 3 2019 동작구 251 408912 211206 18
18 4 2019 금천구 153 251820 122866 13
17 5 2019 구로구 226 439371 219769 17
safe_return_use
10 12602
12 10145
6 20776
20 18766
18 17946
17 15233Loading .shp
seoul.map <- maptools::readShapePoly("./TL_SCCO_SIG_W_SHP/TL_SCCO_SIG_W.shp") # Call .shp file
seoul.nb <- poly2nb(seoul.map) # Builds a neighbours list based on regions with contiguous boundaries
seoul.listw <- nb2listw(seoul.nb) # Supplements a neighbours list with spatial weights for the chosen coding scheme
seoul.mat <- nb2mat(seoul.nb) # Generates a weights matrix for a neighbours list with spatial weights for the chosen coding scheme
# Object of class "nb"
Moran’s I Test
Spatial weight matrix을 이용하여 공간 상관관계에 대한 Moran’s Test를 실시하였다.
\(H_{0}\) : 공간적으로 독립이다. 즉, 공간상관관계가 존재하지 않는다. \(H_{1}\) : 공간상관관계가 존재한다.
### moran.test(observation, listw : a listw object created for example by nb2listw)
moran.test(dat.2019$rape,seoul.listw)
Moran I test under randomisation
data: dat.2019$rape
weights: seoul.listw
Moran I statistic standard deviate = 2.1969, p-value = 0.01401
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.22560008 -0.04166667 0.01480007 ➤ \(p\)값이 0.01401로 유의수준 5%하에서 \(p\)값이 0.05보다 작기 때문에 귀무가설을 기각한다. 즉, 강간강제추행에 대하여 공간상관관계가 존재한다.
Expected Case
흔하게 사용하는 예상되는 수 (\(E_{i}\))는 \[\begin{align*} E_{i} = n_{i}\frac{\sum y_{i}}{\sum n_{i}}, \end{align*}\]
- \(n_{i}\) : \(i\) 번째 지역의 인구수
(참고 : Bayesian Disease Mapping: Hierarchical Modeling in Spatial Epidemiology p.85)
Adjacency Matrix
INLA에서는 근접행렬을 이용하여 그래프로 나타낼 수 있다.
H <- inla.read.graph(filename="Seoul.graph")
image(inla.graph2matrix(H),xlab="",ylab="")
BYM Model
BYM model
\[\begin{align*} \eta_{i} &= \log{\rho_{i}} = b_{0} + u_{i} + v_{i},\\ \hat{\rho_{i}} &= \frac{y_{i}}{E_{i}} = \frac{y_{i}/n_{i}}{\sum y_{i}/ \sum n_{i}}. \end{align*}\]
dat.2019$ID.area <- 1:25 # The Identifiers For The Boroughs
formula <- rape ~ 1 + f(ID.area, model="bym", graph = seoul.adj)
bym <- inla(formula,family = "poisson", data = dat.2019, E = E, # E = E or offset = log(E)
control.predictor=list(compute=TRUE), # Compute the marginals of the model parameters
control.compute = list(dic = TRUE,
waic = TRUE, cpo = TRUE)) # Compute some model choice criteria
summary(bym)
Call:
c("inla(formula = formula, family = \"poisson\", data =
dat.2019, ", " E = E, control.compute = list(dic = TRUE, waic
= TRUE, cpo = TRUE), ", " control.predictor = list(compute =
TRUE))")
Time used:
Pre = 0.636, Running = 0.701, Post = 0.0794, Total = 1.42
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.066 0.105 -0.274 -0.066 0.142 -0.066 0
Random effects:
Name Model
ID.area BYM model
Model hyperparameters:
mean sd 0.025quant
Precision for ID.area (iid component) 3.89 1.11 2.10
Precision for ID.area (spatial component) 1812.27 1815.88 117.06
0.5quant 0.975quant mode
Precision for ID.area (iid component) 3.77 6.42 3.54
Precision for ID.area (spatial component) 1271.56 6628.91 315.90
Expected number of effective parameters(stdev): 24.55(0.125)
Number of equivalent replicates : 1.02
Deviance Information Criterion (DIC) ...............: 231.35
Deviance Information Criterion (DIC, saturated) ....: 49.63
Effective number of parameters .....................: 24.58
Watanabe-Akaike information criterion (WAIC) ...: 224.22
Effective number of parameters .................: 12.57
Marginal log-Likelihood: -155.91
CPO and PIT are computed
Posterior marginals for the linear predictor and
the fitted values are computed# Mean of exp(Fixed effect)
b0 <- bym$marginals.fixed # Extract the marginal posterior distribution for each element of the fixed effects
exp.mean <- lapply(b0,function(x) inla.emarginal(exp,x)) # The exponential transformation and calculate the posterior mean for each of them
exp.mean
$`(Intercept)`
[1] 0.9413661# Confidence interval
exp.ci <- lapply(b0, function(x) inla.qmarginal(c(0.025,0.975),
inla.tmarginal(exp,x))) # inla.tmarginal : transform marginal distribution
exp.ci
$`(Intercept)`
[1] 0.7608768 1.1506961Mapping for spatial effect
Random Effect
# Random Effect for spatial structure (ui+vi)
csi <- bym$marginals.random$ID.area[1:25] # Extract the marginal posterior distribution for each element of the random effects
zeta <- lapply(csi,function(x) inla.emarginal(exp,x)) # The exponential transformation and calculate the posterior mean for each of them
# Define the cutoff for zeta
zeta.cutoff <- c(0.4, 0.8, 1.0, 1.6, 2.0, 2.4, 2.8)
# Transform in categorical variable
cat.zeta <- cut(unlist(zeta),breaks=zeta.cutoff,
include.lowest=TRUE)
➤ BYM Model은 \(\xi_{i}=u_{i}+v_{i}\) 와 \(u_{i}\)를 모수화기 때문에 marginals.random$ID.area의 총 2K개가 존재한다.
➤ 우리가 관심있는 것은 \(\xi_{i}=u_{i}+v_{i}\)이므로 marginals.random$ID.area[1:25]를 추출한다.
Exceedance probabilities
단순히 공간효과만 비교하는 것보다 초과 확률(Exceedance probabilities)을 비교하는 것도 용이하다. 왜냐하면 초과 확률은 불확실성(Uncertainty)을 고려하기 때문이다. 또한, 단순히 점추정이 1이 넘는 지역이라도 불확실성을 고려한 posterior probability를 보면 cutoff value를 넘지 못할 수 있으며, 그 지역이 위험한 지역이라 선언하기 충분한 신뢰성이 없을 수 있다.
➤ \(P(\exp(\xi_{i}>1)\vert Y) = P(\xi_{i}>0\vert Y)\).
maps.cat.zeta <- data.frame(ESRI_PK=dat.2019$ESRI_PK, cat.zeta=cat.zeta, cat.prob=cat.prob)
#Add the categorized zeta to the spatial polygon
data.boroughs <- attr(seoul.map, "data")
attr(seoul.map, "data") <- merge(data.boroughs, maps.cat.zeta,
by="ESRI_PK")
# For adding name of Seoul
lbls <- as.character(seoul.map$SIG_ENG_NM) # 자치구명
spl <- list('sp.text', coordinates(seoul.map), lbls, cex=.7)
# Map zeta
spplot(obj=seoul.map, zcol= "cat.zeta", sp.layout = spl,
col.regions=brewer.pal(6,"Blues"), asp=1)
➤ \(\exp(\xi_{i})\)가 1 보다 크면 근사적으로 \(\rho_{i} > 1\)이므로 \(i\)번째 지역은 전체 지역보다 높은 위험 (risk) 을 가진다.
- 한강을 기준으로 강남 지역이 강북 지역보다 높은 위험을 가지는 지역이 많다.
- 강북 지역의 광진, 용산, 마포, 중구, 종로를 제외한 지역은 전체지역보다 낮은 위험을 가진다.
➤ 가까운 지역들끼리 비슷한 색깔을 가지고 있기 때문에 공간 의존성이 있다고 생각한다.
# Map prob
spplot(obj=seoul.map, zcol= "cat.prob",sp.layout = spl,
col.regions=brewer.pal(5,"Blues"), asp=1)
Estimated Relative Risk
공간효과를 고려하여 추정된 상대적 위험(Relative Risk)는 summary.fitted.values을 이용하여 확인할 수 있다.
# Estimated Relative risk
est.rr <- bym$summary.fitted.values$mean
rr.cutoff <- c(0.4, 0.8, 1.0, 1.6, 2.0, 2.4)
cat.rr <- cut(unlist(est.rr),breaks = rr.cutoff,
include.lowest = TRUE)
# Posterior probability p(Relative risk > a1|y)
mar.rr <- bym$marginals.fitted.values
a1 <- 1
prob.rr <- lapply(mar.rr, function(x) {1 - inla.pmarginal(a1, x)})
cat.rr.prob <- cut(unlist(prob.rr),breaks = prob.cutoff,
include.lowest = TRUE)
maps.cat.rr <- data.frame(ESRI_PK=dat.2019$ESRI_PK, cat.rr=cat.rr, cat.rr.prob=cat.rr.prob)
#Add the categorized zeta to the spatial polygon
data.boroughs <- attr(seoul.map, "data")
attr(seoul.map, "data") <- merge(data.boroughs, maps.cat.rr,
by="ESRI_PK")
# Map
spplot(obj=seoul.map, zcol= "cat.rr",sp.layout = spl,
col.regions=brewer.pal(5,"Blues"), asp=1)
➤ \(\rho_{i}\)가 1 보다 크면, 공간효과를 고려했을 때 \(i\)번째 지역은 전체 지역보다 높은 위험 (risk) 을 가진다.
- 한강을 기준으로 강남 지역이 강북 지역보다 높은 위험을 가지는 지역이 많다.
- 강북 지역의 광진, 용산, 마포, 중구, 종로를 제외한 지역은 전체지역보다 낮은 위험을 가진다.
# Map prob
spplot(obj=seoul.map, zcol= "cat.rr.prob",sp.layout = spl,
col.regions=brewer.pal(5,"Blues"), asp=1)
➤ 0.8 이상인 지역을 위험이 높은 hot-spot 지역, 0.2 이하인 지역을 위험이 낮은 cool-spot 지역이라 한다.
BYM model (Change prior)
BYM Model의 모수 \(\log{\tau_{u}}\)와 \(\log{\tau_{v}}\) 에 대한 기본값은 logGamma(1,0.0005)이다. 모수에 대한 prior을 변경하는 방법은 다음과 같다.
formula1 <- rape ~ 1 + f(ID.area, model="bym", graph = seoul.adj,
hyper = list(prec.unstruct = # Prior for the log of the unstructured effect
list(prior="loggamma",param=c(1,0.01)),
prec.spatial = # Prior for the log of the structured effect
list(prior ="loggamma",param=c(1,0.01))))
bym1 <- inla(formula1,family = "poisson", data = dat.2019, E = E, # E = E or offset = log(E)
control.predictor=list(compute=TRUE), # Compute the marginals of the model parameters
control.compute = list(dic = TRUE,
waic = TRUE, cpo = TRUE)) # Compute some model choice criteria
summary(bym1)
Call:
c("inla(formula = formula1, family = \"poisson\", data =
dat.2019, ", " E = E, control.compute = list(dic = TRUE, waic
= TRUE, cpo = TRUE), ", " control.predictor = list(compute =
TRUE))")
Time used:
Pre = 0.373, Running = 0.675, Post = 0.0843, Total = 1.13
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.066 0.096 -0.258 -0.066 0.125 -0.066 0
Random effects:
Name Model
ID.area BYM model
Model hyperparameters:
mean sd 0.025quant
Precision for ID.area (iid component) 4.14 1.22 2.21
Precision for ID.area (spatial component) 64.15 74.03 1.79
0.5quant 0.975quant mode
Precision for ID.area (iid component) 4.00 6.95 3.71
Precision for ID.area (spatial component) 39.51 264.90 3.13
Expected number of effective parameters(stdev): 24.53(0.131)
Number of equivalent replicates : 1.02
Deviance Information Criterion (DIC) ...............: 231.34
Deviance Information Criterion (DIC, saturated) ....: 49.62
Effective number of parameters .....................: 24.57
Watanabe-Akaike information criterion (WAIC) ...: 224.21
Effective number of parameters .................: 12.57
Marginal log-Likelihood: -152.58
CPO and PIT are computed
Posterior marginals for the linear predictor and
the fitted values are computed