Analisis Regresi Logistik Spasial pada Kejadian Tsunami
2026-10-19
Persiapan dan analisis
Dokumen ini menjalankan analisis regresi logistik biasa dan regresi logistik dengan komponen spasial pada data kejadian tsunami.
# Membersihkan objek lama
rm(list = ls())
# Memanggil library
library(readr)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.5.3library(geosphere)## Warning: package 'geosphere' was built under R version 4.5.3# 1. Membaca data
dataregresilogistikspasial <- read_delim(
"C:/Users/Hype G12/Downloads/dataregresilogistikspasial.csv",
delim = ";",
show_col_types = FALSE
)
# Mengubah data menjadi data frame
data <- as.data.frame(dataregresilogistikspasial)
data$id <- 1:nrow(data)
# 2. Memeriksa data
str(data)## 'data.frame': 782 obs. of 14 variables:
## $ magnitude: num 7 6.9 7 7.3 6.6 7 6.8 6.7 6.8 7.6 ...
## $ cdi : num 8 4 3 5 0 4 1 7 8 9 ...
## $ mmi : num 7 4 3 5 2 3 3 6 7 8 ...
## $ sig : num 768 735 755 833 670 ...
## $ nst : num 117 99 147 149 131 142 136 145 175 271 ...
## $ dmin : num 0.509 2.229 3.125 1.865 4.998 ...
## $ gap : num 17 34 18 21 27 26 22 37 92 69 ...
## $ depth : num 14 25 579 37 624 ...
## $ latitude : num -98 -49.6 -200.5 -192.9 -255.9 ...
## $ longitude: num 160 101 -178 -172 178 ...
## $ Year : num 2022 2022 2022 2022 2022 ...
## $ Month : num 11 11 11 11 11 11 11 10 9 9 ...
## $ tsunami : num 1 0 1 1 1 1 1 1 1 1 ...
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...summary(data)## magnitude cdi mmi sig
## Min. :6.500 Min. :0.000 Min. :1.000 Min. : 650.0
## 1st Qu.:6.600 1st Qu.:0.000 1st Qu.:5.000 1st Qu.: 691.0
## Median :6.800 Median :5.000 Median :6.000 Median : 754.0
## Mean :6.941 Mean :4.334 Mean :5.964 Mean : 870.1
## 3rd Qu.:7.100 3rd Qu.:7.000 3rd Qu.:7.000 3rd Qu.: 909.8
## Max. :9.100 Max. :9.000 Max. :9.000 Max. :2910.0
## nst dmin gap depth
## Min. : 0.0 Min. : 0.00 Min. : 0.00 Min. : 2.70
## 1st Qu.: 0.0 1st Qu.: 0.00 1st Qu.: 14.62 1st Qu.: 14.00
## Median :140.0 Median : 0.00 Median : 20.00 Median : 26.30
## Mean :230.3 Mean : 1.53 Mean : 25.04 Mean : 75.88
## 3rd Qu.:445.0 3rd Qu.: 1.87 3rd Qu.: 30.00 3rd Qu.: 49.75
## Max. :934.0 Max. :161.61 Max. :239.00 Max. :670.81
## latitude longitude Year Month
## Min. :-618.484 Min. :-997.88 Min. :2001 Min. : 1.000
## 1st Qu.: -43.339 1st Qu.: -88.17 1st Qu.:2007 1st Qu.: 3.250
## Median : -3.258 Median : 123.10 Median :2013 Median : 7.000
## Mean : 7.491 Mean : 11.43 Mean :2012 Mean : 6.564
## 3rd Qu.: 38.401 3rd Qu.: 153.19 3rd Qu.:2017 3rd Qu.:10.000
## Max. : 716.312 Max. : 986.03 Max. :2022 Max. :12.000
## tsunami id
## Min. :0.0000 Min. : 1.0
## 1st Qu.:0.0000 1st Qu.:196.2
## Median :0.0000 Median :391.5
## Mean :0.3887 Mean :391.5
## 3rd Qu.:1.0000 3rd Qu.:586.8
## Max. :1.0000 Max. :782.0names(data)## [1] "magnitude" "cdi" "mmi" "sig" "nst" "dmin"
## [7] "gap" "depth" "latitude" "longitude" "Year" "Month"
## [13] "tsunami" "id"table(data$tsunami, useNA = "ifany")##
## 0 1
## 478 304nrow(data)## [1] 782# Memastikan variabel penelitian bertipe numerik
variabel <- c(
"tsunami", "sig", "nst", "gap",
"latitude", "longitude"
)
data[variabel] <- lapply(data[variabel], function(x) {
as.numeric(as.character(x))
})
# Memeriksa data hilang
colSums(is.na(data[variabel]))## tsunami sig nst gap latitude longitude
## 0 0 0 0 0 0# 3. Menyiapkan data untuk analisis
# Memilih variabel yang digunakan
data_analisis <- data[, variabel]
# Menghapus observasi yang memiliki data hilang
data_analisis <- data_analisis[
complete.cases(data_analisis), ]
# Memastikan koordinat berada dalam rentang yang valid
data_analisis <- subset(
data_analisis,
longitude >= -180 & longitude <= 180 &
latitude >= -90 & latitude <= 90
)
# Memastikan respons bernilai 0 atau 1
data_analisis <- subset(
data_analisis,
tsunami %in% c(0, 1)
)
# Memeriksa jumlah observasi dan kategori respons
nrow(data_analisis)## [1] 494table(data_analisis$tsunami)##
## 0 1
## 402 92# Memastikan kedua kategori tsunami tersedia
if (length(unique(data_analisis$tsunami)) != 2) {
stop("Variabel tsunami harus memiliki nilai 0 dan 1.")
}
# 4. Matriks jarak geografis (meter)
coords <- cbind(
data_analisis$longitude,
data_analisis$latitude
)
dist_matrix <- geosphere::distm(
coords,
fun = geosphere::distHaversine
)
# 5. Matriks bobot spasial jarak invers
diag(dist_matrix) <- Inf
# Menangani jarak nol akibat koordinat kembar
dist_matrix[dist_matrix == 0] <- 1
spatial_weights <- 1 / dist_matrix
diag(spatial_weights) <- 0
# Normalisasi baris
spatial_weights <- spatial_weights / rowSums(spatial_weights)
# Memeriksa matriks bobot
dim(dist_matrix)## [1] 494 494dim(spatial_weights)## [1] 494 494head(spatial_weights[, 1:min(5, ncol(spatial_weights))])## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.0000000000 0.0008369237 0.0016036138 0.0043146597 0.0020078905
## [2,] 0.0005232462 0.0000000000 0.0006246113 0.0006486494 0.0007970781
## [3,] 0.0001736057 0.0001081568 0.0000000000 0.0001483843 0.0001969131
## [4,] 0.0029385544 0.0007066055 0.0009334925 0.0000000000 0.0019148192
## [5,] 0.0009811543 0.0006229855 0.0008888073 0.0013738455 0.0000000000
## [6,] 0.0011780908 0.0008709824 0.0009400180 0.0018441124 0.0055298596# 6. Menghitung komponen spasial
data_analisis$spatial_effect <- as.vector(
spatial_weights %*% data_analisis$tsunami
)
head(data_analisis$spatial_effect)## [1] 0.21275998 0.12015702 0.02876439 0.47172834 0.15868024 0.17303374# 7. Regresi logistik biasa
model_logistik <- glm(
tsunami ~ sig + nst + gap,
data = data_analisis,
family = binomial(link = "logit")
)
summary(model_logistik)##
## Call:
## glm(formula = tsunami ~ sig + nst + gap, family = binomial(link = "logit"),
## data = data_analisis)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7057340 0.6934535 2.460 0.0139 *
## sig -0.0019667 0.0008449 -2.328 0.0199 *
## nst -0.0069910 0.0007345 -9.518 <2e-16 ***
## gap 0.0070432 0.0052386 1.344 0.1788
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 474.95 on 493 degrees of freedom
## Residual deviance: 327.70 on 490 degrees of freedom
## AIC: 335.7
##
## Number of Fisher Scoring iterations: 6# 8. Regresi logistik spasial
model_spatial <- glm(
tsunami ~ sig + nst + gap + spatial_effect,
data = data_analisis,
family = binomial(link = "logit")
)
summary(model_spatial)##
## Call:
## glm(formula = tsunami ~ sig + nst + gap + spatial_effect, family = binomial(link = "logit"),
## data = data_analisis)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3912094 0.7845851 -0.499 0.6180
## sig -0.0015082 0.0008912 -1.692 0.0906 .
## nst -0.0056298 0.0007664 -7.346 2.04e-13 ***
## gap 0.0056753 0.0067740 0.838 0.4021
## spatial_effect 7.1254771 1.2475141 5.712 1.12e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 474.95 on 493 degrees of freedom
## Residual deviance: 286.30 on 489 degrees of freedom
## AIC: 296.3
##
## Number of Fisher Scoring iterations: 6# 9. Odds ratio
exp(coef(model_spatial))## (Intercept) sig nst gap spatial_effect
## 0.6762385 0.9984929 0.9943860 1.0056915 1243.2411183# Membuat tabel hasil pendugaan parameter
koef <- summary(model_spatial)$coefficients
tabel_hasil <- data.frame(
Parameter = rownames(koef),
Koefisien = koef[, "Estimate"],
Std_Error = koef[, "Std. Error"],
Z_value = koef[, "z value"],
P_value = koef[, "Pr(>|z|)"],
Odds_Ratio = exp(koef[, "Estimate"]),
row.names = NULL
)
print(tabel_hasil)## Parameter Koefisien Std_Error Z_value P_value Odds_Ratio
## 1 (Intercept) -0.391209430 0.7845851111 -0.4986195 6.180475e-01 0.6762385
## 2 sig -0.001508244 0.0008911853 -1.6924025 9.056925e-02 0.9984929
## 3 nst -0.005629811 0.0007663650 -7.3461216 2.040407e-13 0.9943860
## 4 gap 0.005675322 0.0067740179 0.8378074 4.021389e-01 1.0056915
## 5 spatial_effect 7.125477054 1.2475141447 5.7117405 1.118265e-08 1243.2411183# 10. Perbandingan model
AIC(model_logistik, model_spatial)## df AIC
## model_logistik 4 335.6999
## model_spatial 5 296.2987anova(
model_logistik,
model_spatial,
test = "Chisq"
)## Analysis of Deviance Table
##
## Model 1: tsunami ~ sig + nst + gap
## Model 2: tsunami ~ sig + nst + gap + spatial_effect
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 490 327.7
## 2 489 286.3 1 41.401 1.24e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 11. Prediksi peluang
data_analisis$predicted_prob <- predict(
model_spatial,
type = "response"
)
# Klasifikasi prediksi dengan ambang 0,5
data_analisis$predicted_class <- ifelse(
data_analisis$predicted_prob >= 0.5, 1, 0
)
# 12. Peta sebaran kejadian tsunami
ggplot(
data_analisis,
aes(x = longitude, y = latitude)
) +
geom_point(
aes(color = factor(tsunami)),
size = 2,
alpha = 0.7
) +
labs(
title = "Sebaran Kejadian Tsunami",
x = "Longitude",
y = "Latitude",
color = "Status Tsunami"
) +
theme_minimal()
# 13. Peta peluang hasil model
ggplot(
data_analisis,
aes(x = longitude, y = latitude)
) +
geom_point(
aes(color = predicted_prob),
size = 2,
alpha = 0.7
) +
scale_color_gradient(
low = "green",
high = "red"
) +
labs(
title = "Peluang Tsunami Berdasarkan Model",
x = "Longitude",
y = "Latitude",
color = "Peluang"
) +
theme_minimal()