dataset
Muthi’ah Iffa Karimah
2025-09-10
library(sp)
library(spData)
library(spDataLarge)
library(tmap)
library(geoR)## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.9-6 (built on 2025-08-29) is now loaded
## --------------------------------------------------------------library(crimedata)
library(gstat)
library(sf)## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUElibrary(gstat)
library(maps)
library(spdep)Contoh Data
1. Data Point Pattern (dataset chorley)
dalam pake spatstat.dta, berisi data lokasi rumah penderita kanker di Chorley-Ribble estuary,
library(spatstat.data)
data("chorley")
summary(chorley)## Length Class Mode
## window 5 owin list
## n 1 -none- numeric
## x 1036 -none- numeric
## y 1036 -none- numeric
## markformat 1 -none- character
## marks 1036 factor numericVisualisasi
# Plot visualisasi sebaran kasus
plot(chorley,
main="Chorley Oral Cancer Cases",
cols=c("red","blue"), pch=c(10,10))## Warning in plot.window(...): "cols" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "cols" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "cols" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "cols" is not a
## graphical parameter## Warning in box(...): "cols" is not a graphical parameter## Warning in title(...): "cols" is not a graphical parameter
Eksplorasi
library(spatstat.data) # data chorley
library(spatstat.geom) # ppp, Window## Loading required package: spatstat.univar## spatstat.univar 3.1-4## spatstat.geom 3.5-0library(spatstat.explore) # Gest, quadrat.test, nndist (biasanya tersedia)## Loading required package: spatstat.random## spatstat.random 3.4-1## Loading required package: nlme## spatstat.explore 3.5-2##
## Attaching package: 'spatstat.explore'## The following object is masked from 'package:gstat':
##
## idw# Hitung intensitas
intensity(chorley)## larynx lung
## 0.1840363 3.1032320# Kernel density
d = density(chorley, sigma = 0.3)
plot(d); contour(d, add=TRUE)
d## real-valued pixel image
## 128 x 128 pixel array (ny, nx)
## enclosing rectangle: [343.45, 366.45] x [410.41, 431.79] km# Quadrat test
q <- quadratcount(chorley, nx=4, ny=4)
plot(chorley); plot(q, add=TRUE)
quadrat.test(chorley, nx=4, ny=4)##
## Chi-squared test of CSR using quadrat counts
##
## data: chorley
## X2 = 953.6, df = 15, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 16 tiles (irregular windows)Hasil uji Quadrat test menunjukkan bahwa p-value < 0.05 yaitu 2.2e-16, artinya dapat disimpulkan sebaran kasus kanker di Chorley tidak sesuai dengan CSR. Ada indikasi bahwa pola kemungkinan besar membentuk klaster.
# K-function
K <- Kest(chorley)
plot(K, legend = FALSE, main = "K-function Chorley")
legend("bottomright",
legend = c(expression(hat(K)[iso](r)),
expression(hat(K)[trans](r)),
expression(hat(K)[bord](r)),
expression(K[pois](r))),
col = c("black","red","green","blue"),
lty = 1, bty="n", cex=0.8)
head(K,10)## Function value object (class 'fv')
## for the function r -> K(r)
## ....................................................................
## Math.label Description
## r r distance argument r
## theo K[pois](r) theoretical Poisson K(r)
## border hat(K)[bord](r) border-corrected estimate of K(r)
## trans hat(K)[trans](r) translation-corrected estimate of K(r)
## iso hat(K)[iso](r) Ripley isotropic correction estimate of K(r)
## ....................................................................
## Default plot formula: .~.x
## where "." stands for 'iso', 'trans', 'border', 'theo'
## Recommended range of argument r: [0, 0.093955]
## Available range of argument r: [0, 0.093955]
## Unit of length: 1 kmGaris estimasi (hitam, merah, hijau) lebih tinggi daripada garis teoritis (biru) artinya hampir pada semua jarak r, jumlah pasangan titik lebih banyak dari yang diharapkan di bawah CSR. Sehingga menunjukkan adanya pola klastering pada titik.
# L-function
L = Lest(chorley)
plot(L, . - r ~ r)
head (L,10)## Function value object (class 'fv')
## for the function r -> L(r)
## ....................................................................
## Math.label Description
## r r distance argument r
## theo L[pois](r) theoretical Poisson L(r)
## border hat(L)[bord](r) border-corrected estimate of L(r)
## trans hat(L)[trans](r) translation-corrected estimate of L(r)
## iso hat(L)[iso](r) Ripley isotropic correction estimate of L(r)
## ....................................................................
## Default plot formula: .~.x
## where "." stands for 'iso', 'trans', 'border', 'theo'
## Recommended range of argument r: [0, 0.093955]
## Available range of argument r: [0, 0.093955]
## Unit of length: 1 kmGaris observasi berada di atas garis teoritis (nol) artinya nilai L9r)-r > 0 jarak antar titik lebih kecil dari distribusi CSR. Dapat disimpulkan adanya klastering.
# CSR envelope
env <- envelope(chorley, Kest, nsim=99)## Generating 99 simulations of CSR ...
## 1, 2, [5:32 remaining] 3,
## [5:50 remaining] 4, [6:06 remaining] 5, [5:58 remaining] 6,
## [6:06 remaining] 7, [6:06 remaining] 8, [6:02 remaining] 9,
## [5:55 remaining] 10, [5:51 remaining] 11, [5:45 remaining] 12,
## [5:40 remaining] 13, [5:33 remaining] 14, [5:31 remaining] 15,
## [5:31 remaining] 16, [5:28 remaining] 17, [5:24 remaining] 18,
## [5:20 remaining] 19, [5:13 remaining] 20, [5:09 remaining] 21,
## [5:05 remaining] 22, [5:01 remaining] 23, [4:59 remaining] 24,
## [4:54 remaining] 25, [4:52 remaining] 26, [4:47 remaining] 27,
## [4:43 remaining] 28, [4:40 remaining] 29, [4:36 remaining] 30,
## [4:33 remaining] 31, [4:28 remaining] 32, [4:24 remaining] 33,
## [4:22 remaining] 34, [4:19 remaining] 35, [4:15 remaining] 36,
## [4:11 remaining] 37, [4:06 remaining] 38, [4:02 remaining] 39,
## [3:58 remaining] 40, [3:53 remaining] 41, [3:50 remaining] 42,
## [3:45 remaining] 43, [3:41 remaining] 44, [3:38 remaining] 45,
## [3:33 remaining] 46, [3:30 remaining] 47, [3:26 remaining] 48,
## [3:21 remaining] 49, [3:17 remaining] 50, [3:12 remaining] 51,
## [3:09 remaining] 52, [3:06 remaining] 53, [3:03 remaining] 54,
## [3:00 remaining] 55, [2:56 remaining] 56, [2:53 remaining] 57,
## [2:49 remaining] 58, [2:46 remaining] 59, [2:42 remaining] 60,
## [2:38 remaining] 61, [2:35 remaining] 62, [2:31 remaining] 63,
## [2:27 remaining] 64, [2:23 remaining] 65, [2:19 remaining] 66,
## [2:16 remaining] 67, [2:12 remaining] 68, [2:08 remaining] 69,
## [2:04 remaining] 70, [2:00 remaining] 71, [1:56 remaining] 72,
## [1:52 remaining] 73, [1:48 remaining] 74, [1:44 remaining] 75,
## [1:40 remaining] 76, [1:36 remaining] 77, [1:32 remaining] 78,
## [1:28 remaining] 79, [1:24 remaining] 80, [1:20 remaining] 81,
## [1:16 remaining] 82, [1:12 remaining] 83, [1:07 remaining] 84,
## [1:03 remaining] 85, [59 sec remaining] 86, [55 sec remaining] 87,
## [50 sec remaining] 88, [46 sec remaining] 89, [42 sec remaining] 90,
## [38 sec remaining] 91, [34 sec remaining] 92, [29 sec remaining] 93,
## [25 sec remaining] 94, [21 sec remaining] 95, [17 sec remaining] 96,
## [13 sec remaining] 97, [8 sec remaining] 98, [4 sec remaining]
## 99.
##
## Done.plot(env, legend=FALSE, main="CSR Envelope (K-function)")
legend("topleft",
legend = c("Observed", "Theoretical CSR", "Simulation Envelope"),
col = c("black", "blue", "grey"),
lty = c(1,1,1), lwd=c(2,1,10),
bty="n", cex=0.8)
Garis observasi (hitam( berada di luar batas atas envelope (abu-abu). Sehingga pola spasial cenderung berbeda dari CSR karena keluar ke atas. Artinya ada klastering pada pola titk data.
# G-Function
G = plot(Gest(chorley, legend = FALSE, main="G-function Chorley"))
legend("bottomright",
legend = c(expression(hat(G)(r)), expression(G[pois](r))),
col = c("black","red"),
lty = 1, bty="n", cex=0.8)
head(G,10)## lty col key label meaning
## km 1 1 km italic(hat(G)[km](r)) Kaplan-Meier estimate of G(r)
## rs 2 2 rs italic(hat(G)[bord](r)) border corrected estimate of G(r)
## han 3 3 han italic(hat(G)[han](r)) Hanisch estimate of G(r)
## theo 4 4 theo italic(G[pois](r)) theoretical Poisson G(r)Garis empirirs (hitam, hijau, merah) berada di atas garis teoritis CSR (biru). Sehingga dapat disimpulkan bahwa jarak tetangga terdekat yang lebih kecil lebih tinggi daripada CSR. Artinya titik-titik cenderung lebih berdekatan (ada indikasi klastering.
# F-Function
F = Fest(chorley)
plot(F, legend=FALSE, main="F-function Chorley")
legend("bottomright",
legend = c(expression(hat(F)(r)), expression(F[pois](r))),
col = c("black","red"),
lty = 1, bty="n", cex=0.8)
head(F,10)## Function value object (class 'fv')
## for the function r -> F(r)
## .....................................................................
## Math.label Description
## r r distance argument r
## theo F[pois](r) theoretical Poisson F(r)
## cs hat(F)[cs](r) Chiu-Stoyan estimate of F(r)
## rs hat(F)[bord](r) border corrected estimate of F(r)
## km hat(F)[km](r) Kaplan-Meier estimate of F(r)
## hazard hat(h)[km](r) Kaplan-Meier estimate of hazard function h(r)
## theohaz h[pois](r) theoretical Poisson hazard h(r)
## .....................................................................
## Default plot formula: .~.x
## where "." stands for 'km', 'rs', 'cs', 'theo'
## Recommended range of argument r: [0, 0.37582]
## Available range of argument r: [0, 0.37582]
## Unit of length: 1 kmGaris empiris berada di bawah garis teoritis CSR. Ruang kosong lebih kecil daripada yang diharapkan CSR. Hal ini menandakan bahwa area di sekitar titik relatif lebih padat (indikasi adanya klastering).
2. Lattice (dataset SIDS/sudden infant death syndrome)
sids dalam paket spData/spdep berisi jumalh kasus suddent Infant Death Syndrome di North Carolina
# panggil data
nc <- st_read(system.file("shape/nc.shp", package="sf"))## Reading layer `nc' from data source
## `C:\Users\MUTHI'AH IFFA\AppData\Local\R\win-library\4.5\sf\shape\nc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27# lihat struktur
head(nc)## Simple feature collection with 6 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
## Geodetic CRS: NAD27
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
## 4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
## 5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
## 6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7
## NWBIR74 BIR79 SID79 NWBIR79 geometry
## 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
## 2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
## 3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
## 4 123 830 2 145 MULTIPOLYGON (((-76.00897 3...
## 5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3...
## 6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3...# plot kasus SIDS (BIR74 = jumlah kelahiran 1974, SID74 = kasus SIDS)
plot(nc["SID74"], main = "SIDS Cases 1974 (North Carolina)")
Eksplorasi
Tahapan eksplorasi:
Statistik deskriptif
Hitung rate (misalnya kasus/populasi).
Ringkasan (mean, median, varian) & peta histogram.
Periksa apakah area kecil punya nilai ekstrem (bias kecil sampel).
Visualisasi spasial
Choropleth map → distribusi spasial per area.
Bandingkan peta raw rate vs smoothed (EB).
Autokorelasi spasial global
Uji Moran’s I / Geary’s C.
Mengetahui apakah ada pola cluster secara umum.
Autokorelasi spasial lokal
LISA (Local Moran’s I), Getis-Ord Gi*.
Untuk identifikasi hotspot / coldspot atau outlier spasial.
Smoothing / Empirical Bayes
Untuk mengurangi variabilitas berlebihan di area dengan populasi kecil.
# mengetahui range kasus kelahiran
summary(nc$SID74)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.00 4.00 6.67 8.25 44.00summary(nc$BIR74)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 248 1077 2180 3300 3936 21588# hitung rate SIDS (per 1000 kelahiran)
nc$rate74 <- nc$SID74 / nc$BIR74 * 1000
summary(nc$rate74)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.084 1.855 2.046 2.604 9.554# histogram
hist(nc$rate74, main="Distribusi SIDS Rate (1974)", xlab="Rate per 1000 births")
Visualisasi Choropleth
tmap_mode("plot")## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`tm_shape(nc) +
tm_fill("rate74", style="quantile", palette="Reds",
title="SIDS Rate per 1000 births (1974)") +
tm_borders()##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
## "brewer.reds"Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.
Autokorelasi spasial global (Moran’s I)
# definisi tetangga antar county
nb <- poly2nb(nc)
lw <- nb2listw(nb, style="W")
# uji Moran
moran.test(nc$rate74, lw)##
## Moran I test under randomisation
##
## data: nc$rate74
## weights: lw
##
## Moran I statistic standard deviate = 3.7801, p-value = 7.839e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.230910449 -0.010101010 0.004065134Autokorelasi Spasial Lokal
# Autokorelasi LISA
localI = localmoran(nc$rate74, lw)
nc$Ii = localI[,1]
nc$Ii_p = localI[,5]# Seluruh nilai Ii dan p-value
data.frame(Ii = nc$Ii, p_value = nc$Ii_p)## Ii p_value
## 1 0.6310747658 0.1234702484
## 2 0.6623095293 0.3610115176
## 3 0.2611270896 0.0507183461
## 4 0.0643501181 0.0644743140
## 5 4.5018068705 0.0004422717
## 6 1.7850000887 0.0697788656
## 7 0.6501843571 0.3696436455
## 8 0.0296939150 0.9343922484
## 9 1.0811055175 0.0587380925
## 10 0.1095145801 0.7926973950
## 11 -0.0234643290 0.5129716213
## 12 -0.3326220771 0.4570193677
## 13 -0.0363897639 0.7138680216
## 14 -0.0313095760 0.8556504994
## 15 -0.0521338851 0.5015590496
## 16 1.3937418097 0.0354175922
## 17 0.1185412318 0.1305859330
## 18 0.3832955931 0.0225450561
## 19 0.5699442233 0.1556374877
## 20 -0.0082998474 0.2703254766
## 21 0.2950206655 0.3563741118
## 22 0.6412056874 0.2482071902
## 23 0.4874934946 0.1610425102
## 24 -0.1301780426 0.3741470593
## 25 0.2834456766 0.3044759012
## 26 0.0597938042 0.6999229225
## 27 -0.0638278003 0.7431443722
## 28 2.4869122095 0.0002721252
## 29 -0.0178686841 0.9442255087
## 30 0.0052179655 0.6255389685
## 31 -0.0148395156 0.2604687851
## 32 0.9925945682 0.1746547062
## 33 0.2114655778 0.2701013517
## 34 0.2078518055 0.0330193055
## 35 0.3758372194 0.5396685170
## 36 -0.5326218023 0.0058498689
## 37 0.1203575847 0.5744049930
## 38 -0.3087885221 0.1403691192
## 39 0.3729374556 0.0871082903
## 40 0.5489628762 0.1037719839
## 41 0.7122310822 0.2550468723
## 42 0.1555643647 0.2600784862
## 43 -0.0172359985 0.9185941785
## 44 -0.4425999185 0.6248795392
## 45 -0.4003776654 0.6757161374
## 46 -0.2018287994 0.1332066786
## 47 0.0364695607 0.7575386668
## 48 0.0820758385 0.6396942374
## 49 0.2523810916 0.2781057004
## 50 0.3591557026 0.3004135274
## 51 0.1998081099 0.2213668242
## 52 0.2064635863 0.4787866685
## 53 -0.0159760613 0.9478854917
## 54 0.0381723813 0.7584387016
## 55 -0.1666824697 0.5629805834
## 56 0.8859224498 0.3240768460
## 57 0.0330786473 0.8082641625
## 58 -1.3422028200 0.0784002777
## 59 0.7902325551 0.3044837965
## 60 -0.0343160508 0.5935882522
## 61 -0.1648462785 0.7647748857
## 62 0.1139452223 0.4956238041
## 63 0.0852411349 0.4223438944
## 64 0.0006819882 0.7979708732
## 65 -0.3963712874 0.3310877751
## 66 0.0070170411 0.9370915178
## 67 -0.0151264916 0.6415541630
## 68 0.0016228526 0.4370359636
## 69 0.3675401957 0.3087617675
## 70 0.0864181263 0.2757108644
## 71 0.0183607438 0.3474830755
## 72 -0.0030449304 0.9175451486
## 73 -0.0709721394 0.9424307311
## 74 0.1682790190 0.3666403265
## 75 -0.1224504017 0.4529449583
## 76 -0.1758507118 0.5670825131
## 77 -0.0643641374 0.4011059514
## 78 0.3483000221 0.5214885013
## 79 -0.1140639082 0.5100308214
## 80 -0.0296494635 0.7448119499
## 81 0.0819495254 0.0228202373
## 82 -0.0688240621 0.1344258010
## 83 -0.0271626429 0.7629446019
## 84 -0.6542836391 0.0441278429
## 85 -0.9209450430 0.7405134066
## 86 0.5288285698 0.4459789331
## 87 0.1109301755 0.8414880827
## 88 -0.0630553439 0.5745284363
## 89 -0.5737995453 0.0006084103
## 90 0.8947221534 0.3193849875
## 91 0.0174120694 0.6859444442
## 92 0.5780156639 0.2138929087
## 93 0.0275094997 0.8692631388
## 94 1.3588227417 0.0090064418
## 95 0.0014228794 0.8796069011
## 96 0.9223250256 0.1629222607
## 97 0.3280465193 0.2394819314
## 98 1.4355864412 0.0544826313
## 99 0.0374857093 0.5091837467
## 100 0.1266135698 0.1613305372# Lihat semua atribut nc, termasuk Ii dan Ii_p
head(nc)## Simple feature collection with 6 features and 17 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
## Geodetic CRS: NAD27
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
## 4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
## 5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
## 6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7
## NWBIR74 BIR79 SID79 NWBIR79 geometry rate74
## 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3... 0.9165903
## 2 10 542 3 12 MULTIPOLYGON (((-81.23989 3... 0.0000000
## 3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3... 1.5683814
## 4 123 830 2 145 MULTIPOLYGON (((-76.00897 3... 1.9685039
## 5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3... 6.3335679
## 6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3... 4.8209366
## Ii Ii_p
## 1 0.63107477 0.1234702484
## 2 0.66230953 0.3610115176
## 3 0.26112709 0.0507183461
## 4 0.06435012 0.0644743140
## 5 4.50180687 0.0004422717
## 6 1.78500009 0.0697788656# Hanya tampilkan wilayah dengan autokorelasi signifikan (misal p < 0.05)
subset(data.frame(Ii = nc$Ii, p_value = nc$Ii_p), p_value < 0.05)## Ii p_value
## 5 4.50180687 0.0004422717
## 16 1.39374181 0.0354175922
## 18 0.38329559 0.0225450561
## 28 2.48691221 0.0002721252
## 34 0.20785181 0.0330193055
## 36 -0.53262180 0.0058498689
## 81 0.08194953 0.0228202373
## 84 -0.65428364 0.0441278429
## 89 -0.57379955 0.0006084103
## 94 1.35882274 0.0090064418# peta LISA (nilai Ii)
tm_shape(nc) +
tm_fill("Ii", style="quantile", palette="Blues",
title="Local Moran's I (1974)") +
tm_borders()## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
## "brewer.blues"
## Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.
# Hitung LISA
localI <- localmoran(nc$rate74, lw)
# Simpan nilai Ii dan p-value
nc$Ii <- localI[,1]
nc$Ii_p <- localI[,5]
# Hitung nilai standar (z-score) dari variabel
rate74_mean <- mean(nc$rate74, na.rm = TRUE)
rate74_sd <- sd(nc$rate74, na.rm = TRUE)
z_value <- (nc$rate74 - rate74_mean) / rate74_sd
# Hitung rata-rata tetangga (spatial lag)
lag_rate74 <- lag.listw(lw, nc$rate74)
lag_z <- (lag_rate74 - rate74_mean) / rate74_sd
# Klasifikasi LISA cluster
nc$LISA_cluster <- "NS" # default: Not Significant
alpha <- 0.05 # level signifikansi
nc$LISA_cluster[nc$Ii_p < alpha & z_value > 0 & lag_z > 0] <- "High-High"
nc$LISA_cluster[nc$Ii_p < alpha & z_value < 0 & lag_z < 0] <- "Low-Low"
nc$LISA_cluster[nc$Ii_p < alpha & z_value > 0 & lag_z < 0] <- "High-Low"
nc$LISA_cluster[nc$Ii_p < alpha & z_value < 0 & lag_z > 0] <- "Low-High"
# Tampilkan hasil (10 baris pertama)
head(data.frame(Wilayah = row.names(nc), Ii = nc$Ii, p_value = nc$Ii_p, Cluster = nc$LISA_cluster), 10)## Wilayah Ii p_value Cluster
## 1 1 0.63107477 0.1234702484 NS
## 2 2 0.66230953 0.3610115176 NS
## 3 3 0.26112709 0.0507183461 NS
## 4 4 0.06435012 0.0644743140 NS
## 5 5 4.50180687 0.0004422717 High-High
## 6 6 1.78500009 0.0697788656 NS
## 7 7 0.65018436 0.3696436455 NS
## 8 8 0.02969391 0.9343922484 NS
## 9 9 1.08110552 0.0587380925 NS
## 10 10 0.10951458 0.7926973950 NS# Tampilkan seluruh hasil dengan nama wilayah
hasil_LISA <- data.frame(
Wilayah = nc$NAME, # kolom nama wilayah
Ii = nc$Ii,
p_value = nc$Ii_p,
Cluster = nc$LISA_cluster
)
# Tampilkan semua baris
hasil_LISA## Wilayah Ii p_value Cluster
## 1 Ashe 0.6310747658 0.1234702484 NS
## 2 Alleghany 0.6623095293 0.3610115176 NS
## 3 Surry 0.2611270896 0.0507183461 NS
## 4 Currituck 0.0643501181 0.0644743140 NS
## 5 Northampton 4.5018068705 0.0004422717 High-High
## 6 Hertford 1.7850000887 0.0697788656 NS
## 7 Camden 0.6501843571 0.3696436455 NS
## 8 Gates 0.0296939150 0.9343922484 NS
## 9 Warren 1.0811055175 0.0587380925 NS
## 10 Stokes 0.1095145801 0.7926973950 NS
## 11 Caswell -0.0234643290 0.5129716213 NS
## 12 Rockingham -0.3326220771 0.4570193677 NS
## 13 Granville -0.0363897639 0.7138680216 NS
## 14 Person -0.0313095760 0.8556504994 NS
## 15 Vance -0.0521338851 0.5015590496 NS
## 16 Halifax 1.3937418097 0.0354175922 High-High
## 17 Pasquotank 0.1185412318 0.1305859330 NS
## 18 Wilkes 0.3832955931 0.0225450561 Low-Low
## 19 Watauga 0.5699442233 0.1556374877 NS
## 20 Perquimans -0.0082998474 0.2703254766 NS
## 21 Chowan 0.2950206655 0.3563741118 NS
## 22 Avery 0.6412056874 0.2482071902 NS
## 23 Yadkin 0.4874934946 0.1610425102 NS
## 24 Franklin -0.1301780426 0.3741470593 NS
## 25 Forsyth 0.2834456766 0.3044759012 NS
## 26 Guilford 0.0597938042 0.6999229225 NS
## 27 Alamance -0.0638278003 0.7431443722 NS
## 28 Bertie 2.4869122095 0.0002721252 High-High
## 29 Orange -0.0178686841 0.9442255087 NS
## 30 Durham 0.0052179655 0.6255389685 NS
## 31 Nash -0.0148395156 0.2604687851 NS
## 32 Mitchell 0.9925945682 0.1746547062 NS
## 33 Edgecombe 0.2114655778 0.2701013517 NS
## 34 Caldwell 0.2078518055 0.0330193055 Low-Low
## 35 Yancey 0.3758372194 0.5396685170 NS
## 36 Martin -0.5326218023 0.0058498689 Low-High
## 37 Wake 0.1203575847 0.5744049930 NS
## 38 Madison -0.3087885221 0.1403691192 NS
## 39 Iredell 0.3729374556 0.0871082903 NS
## 40 Davie 0.5489628762 0.1037719839 NS
## 41 Alexander 0.7122310822 0.2550468723 NS
## 42 Davidson 0.1555643647 0.2600784862 NS
## 43 Burke -0.0172359985 0.9185941785 NS
## 44 Washington -0.4425999185 0.6248795392 NS
## 45 Tyrrell -0.4003776654 0.6757161374 NS
## 46 McDowell -0.2018287994 0.1332066786 NS
## 47 Randolph 0.0364695607 0.7575386668 NS
## 48 Chatham 0.0820758385 0.6396942374 NS
## 49 Wilson 0.2523810916 0.2781057004 NS
## 50 Rowan 0.3591557026 0.3004135274 NS
## 51 Pitt 0.1998081099 0.2213668242 NS
## 52 Catawba 0.2064635863 0.4787866685 NS
## 53 Buncombe -0.0159760613 0.9478854917 NS
## 54 Johnston 0.0381723813 0.7584387016 NS
## 55 Haywood -0.1666824697 0.5629805834 NS
## 56 Dare 0.8859224498 0.3240768460 NS
## 57 Beaufort 0.0330786473 0.8082641625 NS
## 58 Swain -1.3422028200 0.0784002777 NS
## 59 Greene 0.7902325551 0.3044837965 NS
## 60 Lee -0.0343160508 0.5935882522 NS
## 61 Rutherford -0.1648462785 0.7647748857 NS
## 62 Wayne 0.1139452223 0.4956238041 NS
## 63 Harnett 0.0852411349 0.4223438944 NS
## 64 Cleveland 0.0006819882 0.7979708732 NS
## 65 Lincoln -0.3963712874 0.3310877751 NS
## 66 Jackson 0.0070170411 0.9370915178 NS
## 67 Moore -0.0151264916 0.6415541630 NS
## 68 Mecklenburg 0.0016228526 0.4370359636 NS
## 69 Cabarrus 0.3675401957 0.3087617675 NS
## 70 Montgomery 0.0864181263 0.2757108644 NS
## 71 Stanly 0.0183607438 0.3474830755 NS
## 72 Henderson -0.0030449304 0.9175451486 NS
## 73 Graham -0.0709721394 0.9424307311 NS
## 74 Lenoir 0.1682790190 0.3666403265 NS
## 75 Transylvania -0.1224504017 0.4529449583 NS
## 76 Gaston -0.1758507118 0.5670825131 NS
## 77 Polk -0.0643641374 0.4011059514 NS
## 78 Macon 0.3483000221 0.5214885013 NS
## 79 Sampson -0.1140639082 0.5100308214 NS
## 80 Pamlico -0.0296494635 0.7448119499 NS
## 81 Cherokee 0.0819495254 0.0228202373 Low-Low
## 82 Cumberland -0.0688240621 0.1344258010 NS
## 83 Jones -0.0271626429 0.7629446019 NS
## 84 Union -0.6542836391 0.0441278429 Low-High
## 85 Anson -0.9209450430 0.7405134066 NS
## 86 Hoke 0.5288285698 0.4459789331 NS
## 87 Hyde 0.1109301755 0.8414880827 NS
## 88 Duplin -0.0630553439 0.5745284363 NS
## 89 Richmond -0.5737995453 0.0006084103 Low-High
## 90 Clay 0.8947221534 0.3193849875 NS
## 91 Craven 0.0174120694 0.6859444442 NS
## 92 Scotland 0.5780156639 0.2138929087 NS
## 93 Onslow 0.0275094997 0.8692631388 NS
## 94 Robeson 1.3588227417 0.0090064418 High-High
## 95 Carteret 0.0014228794 0.8796069011 NS
## 96 Bladen 0.9223250256 0.1629222607 NS
## 97 Pender 0.3280465193 0.2394819314 NS
## 98 Columbus 1.4355864412 0.0544826313 NS
## 99 New Hanover 0.0374857093 0.5091837467 NS
## 100 Brunswick 0.1266135698 0.1613305372 NS# Buat tabel jumlah tiap kategori cluster
table(nc$LISA_cluster)##
## High-High Low-High Low-Low NS
## 4 3 3 90# Konversi hasil table ke data frame
data.frame(table(nc$LISA_cluster))## Var1 Freq
## 1 High-High 4
## 2 Low-High 3
## 3 Low-Low 3
## 4 NS 90- HH (High-High) = Wilayah dengan angka kematian bayi
tinggi dan tetangganya juga tinggi
klaster risiko tertinggi : kemungkinan ada faktor lingkungan, sosial, atau pelayanan kesehatan yang buruk (kemiskinan, faskes tidak memadai, akses faskes sulit, dll)
menjadi area prioritas intervensi
- LL (Low-Low) = Wilayah dengan angka kematian bayi
rendah dan tetangganya juga rendah
klaster sehat : kemungkinan wilayah ini tersedia fasilitas kesehatan yang memadai, tingkat pendidikan ibu tinggi, atau tingkat ekonomi yang lebih baik.
menjadi contoh untuk wilayah lain
- HL (High–Low) = Wilayah dengan angka kematian bayi
tinggi, tetapi tetanggannya rendah
wilayah anomali (outlier) : wilayah ini tidak sesuai dengan asumsi umum yaitu wilayah yang berdekatan biasanya memiliki kondisi yang cenderung sama. Kemungkinan ada faktor lokal tertentu yang membuat angka kematian bayi tinggi tetapi tetangganya rendah.
bisa menjadi penelitian lebih lanjut atas anomali yang terjadi
- LH (Low-High) = Wilayah dengan angka kematian bayi
rendah, tetapi tetangganya tinggi
wilayah anomali positif : kemungkinan ada faktor protektif di wilayah tersebut (program kesehatan lebih baik, fasilitas medis unggul, dll)
bisa menjadi model/best practice
- NS (Not Significant) = tidak ada bukti signifikan
atas autokorelasi spasial
- angka kematian di wilayah ini kemungkinan dipengaruhi oleh faktor acak/lokal yang bukan pola regional
Hotspot Getis-Ord Gi
Gi* mengukur hotpots dan coldspots, nilainya berupa z-score:
Gi* positif = z-score tinggi –> wilayah dengan nilai tinggi dikelilingi nilia tinggi ( hotpots)
Gi* negatif = z-score rendah –> wilayah dengan nilai rendah dikelilingi nilai rendah (coldspots)
Gi* mendekati 0 = netral –> tidak ada pola spasial yang jelas
GI* berbeda dengan Local Moran’s I, nilai Gi* berfokus pada konsentrasi tinggi (hotspot) dan konsentrasi rendah (coldspot)
# Deteksi Hotspot (Getis-Ord Gi*)
gi = localG(nc$rate74, lw)
nc$Gi = as.numeric(gi)
tm_shape(nc) +
tm_fill("Gi", style="quantile", palette="PuOr",
title="Getis-Ord Gi* (Hotspot)") +
tm_borders()## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "PuOr" is named
## "brewer.pu_or"
## Multiple palettes called "pu_or" found: "brewer.pu_or", "matplotlib.pu_or". The first one, "brewer.pu_or", is returned.
##
## [scale] tm_fill:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)
- Orange/merah (Gi* tinggi, signifikan positif)
hotspot angka kematian bayi, daerah klaster dengan risiko tinggi.
perlu ada prioritas intervensi kesehatan masyarakat dan program kesehatan ibu dan anak di wilayah-wilayah ini.
- Biru/ungu (Gi* rendah, signifikan negatif)
coldspot angka kemarian bayi, daerah-daerah ini memiliki tingkat kematian bayi yang rendah serta wilayah tetangganya pun memiliki tingkat kematian yang rendah.
bisa jadi area dengan praktik kesehatan baik, layanan medis bagus, atau kondisi sosial-ekonomi yang lebih baik
- Putih/abu (Gi* netral)
- tidak ada konsentrasi yang signifikan, wilayah ini tidak membentuk hotspot/coldspot.
Model Conditional Autoregressive (CAR) dan Simultanoeus Autoregressive
SAR (Simultaneous Autoregressive Model) = model regresi spasial dengan dependensi di error atau lag
CAR (Conditional Autoregressive Model) = model regresi spasial yang biasa didpakai dalam hierarchical bayesian atau distribusi spasial
# Buat neighbors & spatial weights
nb = poly2nb(nc) # tetangga berdasarkan polygon
lw = nb2listw(nb, style = "W") # bobot normalisasi barisSAR
# install jika belum ada
# install.packages("spatialreg")
# lalu panggil library
library(spatialreg)## Loading required package: Matrix##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption# SAR Lag model
sar_lag = lagsarlm(SID74 ~ BIR74, data = nc, listw = lw)
summary(sar_lag)##
## Call:lagsarlm(formula = SID74 ~ BIR74, data = nc, listw = lw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.36299 -1.89288 -0.55397 1.64313 15.42491
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.146101 0.734718 0.1989 0.8424
## BIR74 0.001702 0.000106 16.0556 <2e-16
##
## Rho: 0.13172, LR test value: 2.3616, p-value: 0.12436
## Asymptotic standard error: 0.087249
## z-value: 1.5097, p-value: 0.13112
## Wald statistic: 2.2792, p-value: 0.13112
##
## Log likelihood: -279.0174 for lag model
## ML residual variance (sigma squared): 15.467, (sigma: 3.9328)
## Number of observations: 100
## Number of parameters estimated: 4
## AIC: 566.03, (AIC for lm: 566.4)
## LM test for residual autocorrelation
## test value: 4.3324, p-value: 0.037394# SAR Error model
sar_err = errorsarlm(SID74 ~ BIR74, data = nc, listw = lw)
summary(sar_err)##
## Call:errorsarlm(formula = SID74 ~ BIR74, data = nc, listw = lw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.11912 -1.65956 -0.60426 1.25865 14.93127
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95278533 0.66302569 1.437 0.1507
## BIR74 0.00174099 0.00010243 16.997 <2e-16
##
## Lambda: 0.33993, LR test value: 7.0396, p-value: 0.0079729
## Asymptotic standard error: 0.1253
## z-value: 2.7129, p-value: 0.0066705
## Wald statistic: 7.3596, p-value: 0.0066705
##
## Log likelihood: -276.6784 for error model
## ML residual variance (sigma squared): 14.421, (sigma: 3.7975)
## Number of observations: 100
## Number of parameters estimated: 4
## AIC: 561.36, (AIC for lm: 566.4)# Model CAR (Conditional Autoregressive)
car_mod <- spautolm(SID74 ~ BIR74, data = nc, listw = lw, family = "CAR")## Warning in spautolm(SID74 ~ BIR74, data = nc, listw = lw, family = "CAR"):
## Non-symmetric spatial weights in CAR modelsummary(car_mod)##
## Call: spautolm(formula = SID74 ~ BIR74, data = nc, listw = lw, family = "CAR")
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.91300 -1.40144 -0.27312 1.29985 13.90157
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.46028995 0.67682700 0.6801 0.4965
## BIR74 0.00178057 0.00010247 17.3764 <2e-16
##
## Lambda: 0.59921 LR test value: 7.1305 p-value: 0.0075783
## Numerical Hessian standard error of lambda: 0.20044
##
## Log likelihood: -276.6329
## ML residual variance (sigma squared): 14.109, (sigma: 3.7561)
## Number of observations: 100
## Number of parameters estimated: 4
## AIC: 561.27AIC(sar_lag, sar_err, car_mod)## df AIC
## sar_lag 4 566.0348
## sar_err 4 561.3568
## car_mod 4 561.26593. Geostatik (dataset jura)
berisi data kandungan logam di Jura, Swiss dengan koordinat spasial.
jura = data(jura)
str(jura.pred)## 'data.frame': 259 obs. of 13 variables:
## $ Xloc : num 2.39 2.54 2.81 4.31 4.38 ...
## $ Yloc : num 3.08 1.97 3.35 1.93 1.08 ...
## $ long : num 6.85 6.85 6.86 6.88 6.88 ...
## $ lat : num 47.1 47.1 47.1 47.1 47.1 ...
## $ Landuse: Factor w/ 4 levels "Forest","Pasture",..: 3 2 2 3 3 3 3 3 3 3 ...
## $ Rock : Factor w/ 5 levels "Argovian","Kimmeridgian",..: 3 2 3 2 5 5 5 1 1 3 ...
## $ Cd : num 1.74 1.33 1.61 2.15 1.56 ...
## $ Co : num 9.32 10 10.6 11.92 16.32 ...
## $ Cr : num 38.3 40.2 47 43.5 38.5 ...
## $ Cu : num 25.72 24.76 8.88 22.7 34.32 ...
## $ Ni : num 21.3 29.7 21.4 29.7 26.2 ...
## $ Pb : num 77.4 77.9 30.8 56.4 66.4 ...
## $ Zn : num 92.6 73.6 64.8 90 88.4 ...head(jura.pred)## Xloc Yloc long lat Landuse Rock Cd Co Cr Cu
## 1 2.386 3.077 6.850413 47.13907 Meadow Sequanian 1.740 9.320 38.32 25.72
## 2 2.544 1.972 6.852674 47.12918 Pasture Kimmeridgian 1.335 10.000 40.20 24.76
## 3 2.807 3.347 6.855886 47.14154 Pasture Sequanian 1.610 10.600 47.00 8.88
## 4 4.308 1.933 6.875799 47.12899 Meadow Kimmeridgian 2.150 11.920 43.52 22.70
## 5 4.383 1.081 6.876925 47.12136 Meadow Quaternary 1.565 16.320 38.52 34.32
## 6 3.244 4.519 6.861415 47.15209 Meadow Quaternary 1.145 3.508 40.40 31.28
## Ni Pb Zn
## 1 21.32 77.36 92.56
## 2 29.72 77.88 73.56
## 3 21.40 30.80 64.80
## 4 29.72 56.40 90.00
## 5 26.20 66.40 88.40
## 6 22.04 72.40 75.20Eksplorasi
summary(jura.pred[, 3:9])## long lat Landuse Rock
## Min. :6.828 Min. :47.12 Forest : 33 Argovian :53
## 1st Qu.:6.849 1st Qu.:47.12 Pasture: 56 Kimmeridgian:85
## Median :6.859 Median :47.13 Meadow :165 Sequanian :63
## Mean :6.858 Mean :47.14 Tillage: 5 Portlandian : 3
## 3rd Qu.:6.867 3rd Qu.:47.15 Quaternary :55
## Max. :6.884 Max. :47.16
## Cd Co Cr
## Min. :0.1350 Min. : 1.552 Min. : 8.72
## 1st Qu.:0.6375 1st Qu.: 6.520 1st Qu.:27.44
## Median :1.0700 Median : 9.760 Median :34.84
## Mean :1.3091 Mean : 9.303 Mean :35.07
## 3rd Qu.:1.7150 3rd Qu.:11.980 3rd Qu.:42.22
## Max. :5.1290 Max. :17.720 Max. :67.60# Histogram
hist(jura.pred$Ni, main="Histogram Ni", xlab="Ni concentration")
Histogram menunjukkan distribusi yang tidak sepenuhnya normal cenderung sedikit menjulur ke kanan.Artinya distribusi menunjukkan adanya daerah hotspot dengan Ni tinggi.
# Boxplot
boxplot(jura.pred[, 3:9], las=2, main="Boxplot logam")
# Plot Lokasi
coordinates(jura.pred) <- ~Xloc + Yloc
plot(jura.pred, pch=20, col="blue", main="Lokasi Sampel")
spplot(jura.pred, "Ni", main="Sebaran Ni")
# Hubungan antara variabel
pairs(jura.pred@data[, 3:9], main="Scatterplot Matrix")
# ambil variabel logam saja
logam = jura.pred[, c("Co","Cr","Cu","Ni","Pb","Zn","Cd")]
# pastikan jadi data.frame biasa
logam <- as.data.frame(logam)
# korelasi
cor(logam)## Co Cr Cu Ni Pb Zn
## Co 1.0000000 0.45343574 0.21668062 0.75073685 0.1887371 0.47429323
## Cr 0.4534357 1.00000000 0.21006986 0.69268265 0.2950831 0.67340647
## Cu 0.2166806 0.21006986 1.00000000 0.22937845 0.7783533 0.57005245
## Ni 0.7507369 0.69268265 0.22937845 1.00000000 0.3081438 0.63466772
## Pb 0.1887371 0.29508307 0.77835330 0.30814382 1.0000000 0.59153289
## Zn 0.4742932 0.67340647 0.57005245 0.63466772 0.5915329 1.00000000
## Cd 0.2533525 0.60941668 0.11985435 0.48737516 0.2223784 0.66920421
## Xloc 0.3065004 -0.05333416 -0.04320191 0.23708348 -0.1392969 0.07378990
## Yloc -0.1609440 0.03471288 -0.10255323 -0.09359487 -0.2101584 -0.02867539
## Cd Xloc Yloc
## Co 0.2533525 0.30650041 -0.16094399
## Cr 0.6094167 -0.05333416 0.03471288
## Cu 0.1198543 -0.04320191 -0.10255323
## Ni 0.4873752 0.23708348 -0.09359487
## Pb 0.2223784 -0.13929695 -0.21015840
## Zn 0.6692042 0.07378990 -0.02867539
## Cd 1.0000000 0.11395876 0.11233764
## Xloc 0.1139588 1.00000000 0.19808972
## Yloc 0.1123376 0.19808972 1.00000000# scatterplot matrix
pairs(logam)
vgm_exp <- variogram(Ni ~ 1, jura.pred)
plot(vgm_exp)
# --- Variogram & fitting ---
vgm_exp <- variogram(Ni ~ 1, jura.pred)
vgm_model <- vgm(psill=100, model="Sph", range=5, nugget=10)
vgm_fit <- fit.variogram(vgm_exp, model=vgm_model)
# --- Grid untuk prediksi ---
bb <- bbox(jura.pred)
x.range <- seq(bb["Xloc",1], bb["Xloc",2], length.out=100)
y.range <- seq(bb["Yloc",1], bb["Yloc",2], length.out=100)
grd <- expand.grid(Xloc = x.range, Yloc = y.range)
coordinates(grd) <- ~Xloc+Yloc
gridded(grd) <- TRUE
proj4string(grd) <- proj4string(jura.pred)
# --- Ordinary Kriging ---
krig_res <- krige(Ni ~ 1, jura.pred, grd, model=vgm_fit)## [using ordinary kriging]library(gridExtra)
# Plot prediksi
p1 <- spplot(krig_res["var1.pred"], main="Ordinary Kriging (Ni)")
# Plot variansi
p2 <- spplot(krig_res["var1.var"], main="Var Kriging (Ni)")
# Tampilkan berdampingan
grid.arrange(p1, p2, ncol=2)
Plot Prediksi Ordinary Kriging (Ni)
Kuning –> Ungu : menunjukkan rentang prediksi konsentrasi 5 - 40
Kuning–> Orange : menunjukkan konsentrasi tinggi 30-40
Biru tua : menunjukkan konsentrasi rendah 5-10
Biru ungu –> Hitam : spot tinggi ke rendah yang menunjukkan variasi spasial cukup kuat
Distribusi Ni cenderung tidak merata, ada hotspot dengan nilai tinggi dan area nilai rendah.
Plot Variansi Prediksi Kriging (Ni)
Warna menunjukkan ketidakpastian prediksi (10-90)
Biru gelap –> variansi kecil: prediksi lebih akurat karena dekat dengan titik sampel asli
Kuning –> variansi besar: prediksi cenderung kurang pasti karena jauh dari titik sampel asli