Spatial Analysis

1 Pendahuluan

1.1 Data accidents

1.2 Data accident_data

## class       : SpatialPointsDataFrame 
## features    : 8 
## extent      : -19.99336, -19.81216, -44.02575, -43.92826  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## variables   : 3
## names       : weekDay,  hour,         type 
## min values  :     mon, 00:30,    collision 
## max values  :     sun, 22:50, running over

1.3 Peta Dunia

1.4 Peta Indonesia dengan berbagai macam proyeksi

1.4.1 Proyeksi conic

1.4.2 Proyeksi bonne

1.4.3 Proyeksi albers

1.4.4 Proyeksi lagrange

2 Model Spasial Asuransi Kendaraan Bermotor

Data yang digunakan adalah data kecelakaan kendaraan bermotor di Brazil. Dalam data ini terdapat fitur-fitur berikut:

##  [1] "CityCode"   "State"      "City"       "PopExpo"    "PopPrem"   
##  [6] "PopClaimRo" "PopClaimCo" "PopClaimFi" "PopClaimOt" "LuxExpo"   
## [11] "LuxPrem"    "LuxClaimRo" "LuxClaimCo" "LuxClaimFi" "LuxClaimOt"
## [16] "HDIcity00"  "CityDens10" "StateAb"

Keterangan :

• CityCode = severity untuk setiap kota
  
• State = nama lengkap state
• City = nama lengkap kota
• PopExpo (LuxExpo)= jumlah periode eksposur untuk polis kendaraan populer (mewah), dalam tahun
• PopPrem (LuxPrem)= Gross premium untuk grup kendaraan populer (mewah)
• PopClaimRo (LuxClaimRo) = total banyak klaim perampokan untuk kendaraan populer (mewah)
• PopClaimCo (LuxClaimCo) = total banyak klaim tabrakan untuk kendaraan populer (mewah)
• PopClaimFi (LuxClaimFi) = total banyak klaim kebakaran untuk kendaraan populer (mewah)
• PopClaimOt (LuxClaimOt) = total banyak klaim lainnya untuk kendaraan populer (mewah)
• HDIcity00 = Human Development Index untuk setiap kota
• CityDens10 = kepadatan populasi untuk setiap kota
• StateAb = nama singkatan state ( 2 huruf)

2.0.1 Penggunaan Paket RgoogleMaps

map.susep<-SpatialPolygons2PolySet(shape)
bb<-bbox(shape)
MyMap<-GetMap.bbox(bb[1,],bb[2,],maptype="satellite",destfile = "myMap.png",GRAYSCALE = FALSE)
plotvar<-shape$HDIcity00
ncls<-5

# membuat warna 
colpal<-brewer.pal(ncls,"Greens")

# membagi nilai plotvar kedalam ncls interval dengan style
classes<-classIntervals(plotvar,ncls,style = "equal")

## Warning in classIntervals(plotvar, ncls, style = "equal"): var has missing
## values, omitted in finding classes

#Memetakan kelas ke warnanya
cols2<- findColours(classes,colpal)

# missing value harus diisi
cols2[is.na(shape$HDIcity00)]<-"red"

# Menampilkan map.susep dengan latar belakang MyMap
PlotPolysOnStaticMap(MyMap,map.susep,col=cols2,lwd=0.15,border=NA,add=FALSE)
legend("topleft",fill=attr(cols2,"palette"),legend=leglabs(round(classes$brks,digits=2)),
       cex=1.0,ncol=1,bg="white",bty="o")

#cols <- rev(gray(seq(0.1, 0.9, length = 5)))
# Human Development Index
#spplot(shape, "HDIcity00", col.regions = cols, cuts = length(cols) - 1)
#spplot(brgeomunicins, "PopClaimFire", col.regions = cols, cuts = length(cols) - 1)
#spplot(brgeomunicins, "PopClaimColl", col.regions = cols, cuts = length(cols) - 1)
#spplot(brgeomunicins, "PopClaimRob", col.regions = cols, cuts = length(cols) - 1)

2.1 Peta Ketetanggan suatu State

Disini kami tampilkan beberapa jenis ketetanggan dari suatu state

2.1.1 Ketetanggan berdasarkan garis perbatasan

Kita sebut dua daerah bertetangga bila kedua daerah tersebut memiliki garis perbatasan yang sama.

pos<-which(shape@data$StateAb=="PR")
prshape<-shape[pos,]
#plot(prshape)
#text(coordinates(prshape),label=prshape@data$City,cex=0.5)

# membuat garis ketetanggaan tanpa bobot
pr.nb<-poly2nb(prshape)

# membuat plot poligon state PR
plot(prshape)

#menambahkan garis hubungan ketetanggaan dalam plot poligon PR 
plot(pr.nb,coordinates(prshape),add=TRUE,col="blue")
title("Keteanggaan berdasarkan garis perbatasan")

Visalisasi di atas terlihat sangat kompleks dan tidak mudah untuk diamati. Sebagai alternatif ketetanggaan kita dapat menggunakan definisi ketetanggan berdasarkan tetangga terdekat berdasarkan jarak Euclidean yang ditampilkan pada bagian berikut.

2.1.2 Ketetangaan berdasarkan K nearest neighborhood (knn)

Disini dua daerah yang memiliki garis perbatasan yang sama belum tentu bertetangga. Kita pilih titik pusat suatu daerah sebagai titik acuan dari setiap daerah. \(K\) tetangga terdekat dari daerah \(i\) adalah \(K\) daerah dengan pusat titik terdekat ke pusat daerah \(i\). Disini kita gunakan jarak Euclidean untuk mengukur jarak dari pusat daerah \(i\) ke pusat daerah \(j\).

coords<-coordinates(prshape)
pr.knn<-knearneigh(coords,k=3,longlat=TRUE)
pr.nbknn<-knn2nb(pr.knn)
plot(prshape,border="grey")
plot(pr.nbknn,coords,add=TRUE, col="blue")
title("K nearest neighbours, k=3")

Terlihat pada gambar di atas dengan menggunakan ketetanggaan berdasarkan 3 tetangga terdekat visualisasi yang kita peroleh menjadi lebih sederhana.

2.1.3 Moran’s Index

Moran’s Index ini digunakan untuk mengukur asosiasi spasial dari suatu daerah. Ukuran empirik dari Moran’s Index ini adalah

\[ I=\frac{1}{\sum_{i\neq j} w_{ij}} \sum_{i\neq j}w_{ij}\left(\frac{y_i-\overline{y}}{s_y} \right)\left(\frac{y_j-\overline{y}}{s_y} \right).\] Asumsikan \(Y_1, Y_2,...,Y_n\) saling bebas dan berdistribusi identik. Kami gunakan algoritma berikut:

1. Hitung Moran’s Index \(I\) untuk data observasi \(y_1,y_2, ...,y_n\), namakan \(I^{(1)}\).
2. untuk \(k=1\) sampai \(M\) lakukan

• permutasikan nilai \(y_i\)
• hitung Moran’s Index untuk setiap permutasi ini akan diperoleh \(I^{(k+1)}\)

Kita gunakan hipotesa: \(H_o\): semua permutasi mempunyai peluang yang sama

dengan menggunakan p-value: \[p-\text{value}= \frac{Banyaknya I^{(j)}>I^{(1)},\ j=1,2,...,M}{M+1}.\] kita tolah \(H_o\) jika \(p-\)value \(<\alpha\) dengan \(\alpha\) suatu tingkat signifikansi yang kita pilih ( misalkan \(\alpha =5\%\))

pr.listw<- nb2listw(pr.nb,style="W")
imoran<-moran.mc(prshape$HDIcity00,pr.listw,nsim=999 )
par(mar=c(4,4,2,2))
hist(imoran$res,xlab='Index',main='',col=gray(.5), border=gray(.7))
arrows(imoran$stat,-2,imoran$stat,10,lwd=2,col=2,leng=.1,code=1)
segments(imoran$stat,3,0.4,120,lty=2)
text(.4,150,paste("Moran's I=", format(imoran$stat,dig=4)))
text(.4,130,paste("p-value=", format(imoran$p.val,dig=4)))

Dari hasil di atas kita peroleh \(p\)-value =\(0.001=0.1\%\), sehingga kita dapat menolak \(H_o\). Dalam kasus ini terdapat asosiasi spasial pada data Human Development Index.

2.2 Spatial Car Accident Insurance Analysis

2.2.1 Analisa Spasial untuk masing-masing kategori kendaraan

Disini kita gunakan model Conditional Auto-regression (CAR) untuk masing-masing kategori kendaraan. Dalam model ini pengaruh antar kategori tidak diperhitungkan. Parameter spasial \(\phi\) dimodelkan sebagai berikut:

\[\phi_i | \phi_{-i}\sim N\left(\rho \frac{\sum_{i\sim j}\phi_j}{n_i},\frac{\sigma^2}{n_i}\right),\ i=1,2,...,n,\] dengan \(i\sim j\) menyatakan semua tetangga dari \(i\), \(\sigma^2\) variansi, dan \(\rho\) kekuatan dari hubungan spasial.

Diberikan data \(Y_i\) yang menyatakan banyaknya klaim kendaraan di daerah \(i\). Kita modelkan \[Y_i\sim Poiss(E_i\theta_i),\ i=1,2,...,n,\] \[\log(\theta_i)=\eta_i+\phi_i+e_i,\] dengan \(E_i\) ekspektasi laju banyaknya klaim pada daerah \(i\), \(\eta_i=X^T_i\beta\), \(\phi_i\) efek acak spasial pada daerah \(i\), dan \(e_i\) gangguan acak.

Perhatikan dalam model ini laju banyaknya klaim dapat lebih besar atau lebih kecil dari ekspektasinya. Jika \(\eta_i+\phi_i+e_i>1\) maka lajunya akan melebihi ekspektasinya. Sedangkan jika \(\eta_i+\phi_i+e_i<1\), lajunya akan lebih rendah dari ekspektasinya.

Kita gunakan Human Development Index (HDI) dan logaritma dari kepadatan penduduk sebagai fiturnya dalam \(\eta_i\), \[\eta_i=\beta_0+\beta_1HDI_i+\beta_2 Kepadatan_i,\ i=1,2,...,n.\]

2.3 Pengolahan data hilang

Dalam data terdapat informasi exposure untuk kendaraan popular dan mewah dan data HDI di beberapa lokasi yang hilang. Kita dapat mengisi data yang hilang ini dengan menggunakan regresi linear sederhana.

2.3.1 Pengisian data hilang Human Development Index

# Fungsi rataan
NbMean<-function(shp,vari){
  # Buat daftar tetangga untuk setiap daerah di shp
  shpnb<-poly2nb(shp)
  
  # Buat matriks ketetanggan dengan style bobot binary
  shpnb.mat<-nb2mat(shpnb,style="B",zero.policy=TRUE)
  
  # cari daerah yang nilainya hilang
  selNA<-which(is.na(shp@data[,vari]))
  
  # indeks yang nilainya hilang
  NAnb<-shpnb.mat[selNA,]
  
  # Pengisian rataan pada data hilang
  shp@data[selNA,vari]<- apply(NAnb,1,FUN=function(x) 
    mean(shp@data[which(x==1),vari],na.rm=TRUE))
  return(shp)

}

# Pegisian data hilang HDI
shape<-NbMean(shape,"HDIcity00")

2.3.2 Pengisian data hilang Exposure kendaraan Popular

# Total klaim untuk kendaraan popular
shape@data$SIN_POP<-rowSums(shape@data[,6:9],na.rm=FALSE)

# Total klaim untuk kendaraan mewah
shape@data$SIN_LUX<-rowSums(shape@data[,12:15],na.rm=FALSE)

# logaritma data kepadatan penduduk setiap daerah
shape@data$POP<- log(shape@data$CityDens10)

# Pengecekan data hilang di exposure kendaraan popular
# posisi data yang tidak hilang
pos<-which(!is.na(shape@data$PopExpo))
#posisi data hilang
pos.ms<-which(is.na(shape@data$PopExpo))

# kita logaritmakan data exposure pada posisi yang tidak hilang
expo<-log(shape@data$PopExpo[pos]+1)

# kepadatan pada posisi data exposure popular yang tidak hilang
pop<-shape@data$POP[pos]

#bangun model regresi sederhana dengan respon expo dan predictor pop
reg<-lm(expo~pop)
# hitung nilai regresi pada lokasi exposure popular yang hilang
pred<-as.vector(predict(reg,data.frame(pop=shape@data$POP[pos.ms])))

# kembalikan ke nilai eksposure aslinya
shape@data$PopExpo[pos.ms]<-sapply((exp(pred)-1),function(x) max(x,0))

2.3.3 Pengisian data hilang Exposure kendaraan Mewah

# Pengecekan data hilang di exposure kendaraan mewah
# posisi data yang tidak hilang
pos<-which(!is.na(shape@data$LuxExpo))
#posisi data hilang
pos.ms<-which(is.na(shape@data$LuxExpo))

# kita logaritmakan data exposure pada posisi yang tidak hilang
expo<-log(shape@data$LuxExpo[pos]+1)

# kepadatan pada posisi data exposure mewah yang tidak hilang
pop<-shape@data$POP[pos]

#bangun model regresi sederhana dengan respon expo dan predictor pop
reg<-lm(expo~pop)
# hitung nilai regresi pada lokasi exposure mewah yang hilang
pred<-as.vector(predict(reg,data.frame(pop=shape@data$POP[pos.ms])))

# kembalikan ke nilai eksposure aslinya
shape@data$LuxExpo[pos.ms]<-sapply((exp(pred)-1),function(x) max(x,0))

2.3.4 Ekspektasi banyaknya kecelakaan berdasarkan tipe kendaraan popular(mewah)

Kita hitung proporsi eksposure kendaraan popular(mewah) untuk setiap daerahnya sebagai berikut: \[ \text{Prop}_i^{pop(mwh)}=\frac{\text{banyaknya eksposure kendaraan popular(mewah) daerah ke-} i}{\text{total eksposure kendaraan popular(mewah)}} ,\ i=1,...,n\]

Sehingga ekspektasi banyaknya kendaraan popular(mewah) mengalami kecelakaan: \[p^{pop(mwh)}_i=\text{Prop}_i^{pop(mwh)} \times \text{Banyaknya klaim mobil popular(mewah) di daerah ke-}i,\ i=1,...,n\] Berikut ini perhitungan dalam R:

# ekspektasi banyaknya kecelakaan kendaraan popular
shape@data$E_POP<-shape@data$PopExpo*sum(shape@data$SIN_POP,na.rm=TRUE)/
  sum(shape@data$PopExpo)

#ekspektasi banyaknya kecelakaan kendaraan mewah
shape@data$E_LUX<-shape@data$LuxExpo*sum(shape@data$SIN_LUX,na.rm=TRUE)/
  sum(shape@data$LuxExpo)

2.4 Spatial Bayesian dengan Integrated Nested Laplace Approximation (INLA)

Kita buat struktur spasial untuk penggunaan INLA sebagai berikut:

shape@data$struct<-rep(1:dim(shape@data)[1])
shape@data$unstruct<-rep(1:dim(shape@data)[1])

nb2INLA("ngbINLA.graph",poly2nb(shape))

f.pop<-SIN_POP~HDIcity00 + POP +f(unstruct,model="iid")+
f(struct,model="besag",graph="ngbINLA.graph")

f.lux<-SIN_LUX~HDIcity00 + POP +f(unstruct,model="iid")+
f(struct,model="besag",graph="ngbINLA.graph")


m.pop<-inla(f.pop,family="poisson", data=shape@data,E=E_POP,
            control.compute = list(dic=TRUE,cpo=TRUE),
            control.predictor=list(compute=TRUE,link=1))
summary(m.pop)

## 
## Call:
## c("inla(formula = f.pop, family = \"poisson\", data = shape@data, ",  "    E = E_POP, control.compute = list(dic = TRUE, cpo = TRUE), ",  "    control.predictor = list(compute = TRUE, link = 1))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          4.7419        260.9263          1.8426        267.5108 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -1.3617 0.4027    -2.1558  -1.3607    -0.5743 -1.3585   0
## HDIcity00   -0.2206 0.5931    -1.3834  -0.2213     0.9452 -0.2228   0
## POP          0.1111 0.0149     0.0821   0.1111     0.1404  0.1110   0
## 
## Random effects:
## Name   Model
##  unstruct   IID model 
## struct   Besags ICAR model 
## 
## Model hyperparameters:
##                          mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for unstruct 27.468 7.907     15.519   26.261      46.23 24.02
## Precision for struct    6.506 1.643      3.794    6.338      10.22  6.02
## 
## Expected number of effective parameters(std dev): 365.48(18.77)
## Number of equivalent replicates : 3.929 
## 
## Deviance Information Criterion (DIC) ...............: 5670.06
## Deviance Information Criterion (DIC, saturated) ....: 1946.29
## Effective number of parameters .....................: 364.02
## 
## Marginal log-Likelihood:  -4373.82 
## CPO and PIT are computed
## 
## Posterior marginals for linear predictor and fitted values computed

m.lux<-inla(f.lux,family="poisson", data=shape@data, E=E_LUX,
            control.compute = list(dic=TRUE,cpo=TRUE),
            control.predictor=list(compute=TRUE,link=1))
summary(m.lux)

## 
## Call:
## c("inla(formula = f.lux, family = \"poisson\", data = shape@data, ",  "    E = E_LUX, control.compute = list(dic = TRUE, cpo = TRUE), ",  "    control.predictor = list(compute = TRUE, link = 1))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          3.9218        205.5388          0.8888        210.3493 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -1.0179 0.7211    -2.4434  -1.0148     0.3892 -1.0086   0
## HDIcity00    0.0780 1.0503    -1.9757   0.0748     2.1476  0.0684   0
## POP          0.0507 0.0258     0.0001   0.0507     0.1015  0.0506   0
## 
## Random effects:
## Name   Model
##  unstruct   IID model 
## struct   Besags ICAR model 
## 
## Model hyperparameters:
##                          mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for unstruct 11.031 3.324      5.999   10.532     18.944 9.609
## Precision for struct    5.239 1.925      2.396    4.934      9.877 4.370
## 
## Expected number of effective parameters(std dev): 221.21(18.67)
## Number of equivalent replicates : 6.492 
## 
## Deviance Information Criterion (DIC) ...............: 3311.25
## Deviance Information Criterion (DIC, saturated) ....: NaN
## Effective number of parameters .....................: 218.17
## 
## Marginal log-Likelihood:  -3117.63 
## CPO and PIT are computed
## 
## Posterior marginals for linear predictor and fitted values computed

2.4.1 Visualisasi Kebergantungan Spasial

# Visualisasi

shape@data$POP_SPAT<-m.pop$summary.random$struct$mean
shape@data$POP_PRED<-m.pop$summary.fitted.values$mean*shape@data$E_POP
shape@data$POP_RR<- m.pop$summary.linear.predictor$mean

shape@data$LUX_SPAT<-m.lux$summary.random$struct$mean
shape@data$LUX_PRED<-m.lux$summary.fitted.values$mean*shape@data$E_LUX
shape@data$LUX_RR<- m.lux$summary.linear.predictor$mean

spplot(shape,c("POP_SPAT", "LUX_SPAT"),layout=c(2,1),
       main="spatial dependence",cuts=5,col.region=grey.colors(50,1,0))

Plot kebergantungan spasial dapat digunakan untuk menganalisa potensi daerah dalam kontribusi klaim.

2.4.2 Visualisasi log relatif risiko Spasial

spplot(shape,c("POP_RR", "LUX_RR"),layout=c(2,1),
       main="log relative risk by category",cuts=5,col.region=grey.colors(50,1,0))

Dengan plot spasial di atas kita dapat memperoleh informasi daerha mana saja yang memiliki risiko yang besar untuk kedua kategori kendaraan.

2.4.3 Plot Prediksi vs Observasi

par(mfrow=c(1,2))

plot(log(shape@data$POP_PRED+1),log(shape@data$SIN_POP+1),
     xlab="log pred kecelakaan kendaraan popular", 
     ylab="log data observasi kecelakaan kendaraan popular", 
     main="log-predicted x log-accidents")

abline(a=0,b=1)

plot(log(shape@data$LUX_PRED+1),log(shape@data$SIN_LUX+1),
     xlab="log pred kecelakaan kendaraan mewah", 
     ylab="log data observasi kecelakaan kendaraan mewah", 
     main="log-predicted x log-accidents")

Dapat dilihat dalam plot di atas, hasil prediksi dapat menjelaskan observasi dengan baik untuk setiap kategori kendaraan.

2.4.4 Analisa spasial memperhitungkan efek spasial bersama

Bila kita memperhitungkan efek spasial bersama antara kedua kategori, model CAR yang digunakan sebagai berikut: \[Y_{ji}\sim Poiss(E_{ji}\theta_{ji}), \ j=\{\text{popular, mewah}\}, \ i=1,...,n,\] \[\log(\theta_{ji})=\eta_{ji}+\gamma(j)\psi_{i}+\phi_{ji}+e_{ij},\] dengan \(\eta_{ji}=\beta_0+\beta_1^{(j)}HDI_i^{(j)}+\beta_2^{(j)}POP_i^{(j)}\), \(\psi_i\) adalah komponen spasial bersama dan \[\gamma(j)=\begin{cases} 1,\text{ jika } j=\text{ popular}\\ \delta, \text{ jika } j=\text{ mewah} \end{cases}\]

k<- 2
n<- dim(shape@data)[1]
Y<-matrix(NA,n,k)
Y[1:n,1]<-shape@data$SIN_POP
Y[1:n,2]<- shape@data$SIN_LUX


share.dat<-list(Y=matrix(NA,nrow=n*2,ncol=2))
share.dat$Y[1:n,1]<- Y[,1]
share.dat$Y[n+(1:n),2]<-Y[,2]

share.dat$E<- c(shape@data$E_POP,shape@data$E_LUX)

Untuk analisa lebih lanjut dapat anda lihat pada materi perkulihan yang saya akan berikan semester depan.