SYNTAX ANALISIS STUDI CROSS-SECTIONAL DAN PEMETAAN RISIKO DEMAM BERDARAH DENGUE (DBD) DI PROVINSI JAWA BARAT TAHUN 2024
Nama: Fatih Zahrani
NPM:
140610230014
Dosen Pengampu: I Gede Nyoman Mindra
Jaya, M.Si., Ph.D.
Mata Kuliah: Epidemiologi
S1 Statistika FMIPA UNPAD
Syntax Analisis
library(sf)## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(openxlsx)## Warning: package 'openxlsx' was built under R version 4.3.3library(viridis)## Warning: package 'viridis' was built under R version 4.3.3## Loading required package: viridisLite## Warning: package 'viridisLite' was built under R version 4.3.3library(spdep)## Warning: package 'spdep' was built under R version 4.3.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.3.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`# -------------------------- MASUKKAN PLOT SPASIAL -----------------------------
Indo <- st_read("D:/Kuliah/Sem 5/Epidem/Dashboard/shapefile/RBI_50K_2023_Jawa Barat.shp")## Reading layer `RBI_50K_2023_Jawa Barat' from data source
## `D:\Kuliah\Sem 5\Epidem\Dashboard\shapefile\RBI_50K_2023_Jawa Barat.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 25 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: 106.3703 ymin: -7.82099 xmax: 108.8468 ymax: -5.806538
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84names(Indo)## [1] "NAMOBJ" "FCODE" "REMARK" "METADATA" "SRS_ID"
## [6] "KDBBPS" "KDCBPS" "KDCPUM" "KDEBPS" "KDEPUM"
## [11] "KDPBPS" "KDPKAB" "KDPPUM" "LUASWH" "TIPADM"
## [16] "WADMKC" "WADMKD" "WADMKK" "WADMPR" "WIADKC"
## [21] "WIADKK" "WIADPR" "WIADKD" "SHAPE_Leng" "SHAPE_Area"
## [26] "geometry"Indo$WADMKK## [1] "Bandung" "Bandung Barat" "Bekasi" "Bogor"
## [5] "Ciamis" "Cianjur" "Cirebon" "Garut"
## [9] "Indramayu" "Karawang" "Kota Bandung" "Kota Banjar"
## [13] "Kota Bekasi" "Kota Bogor" "Kota Cimahi" "Kota Cirebon"
## [17] "Kota Depok" "Kota Sukabumi" "Kota Tasikmalaya" "Kuningan"
## [21] "Majalengka" "Pangandaran" "Purwakarta" "Subang"
## [25] "Sukabumi" "Sumedang" "Tasikmalaya"Indo <- st_zm(Indo, drop = TRUE, what = "ZM")
ggplot(Indo) +
geom_sf(fill = "lightblue", color = "black") +
labs(title = "Peta Kabupaten/Kota di Provinsi Jawa Barat") +
theme_minimal()
# ------------------------------- MASUKKAN DATA --------------------------------
file_path <- "D:/Kuliah/Sem 5/Epidem/Dashboard/data/DATA EPIDEMIOLOGI.xlsx"
data <- read.xlsx(file_path)
head(data)## no kab_kota JK JP KDBD
## 1 1 Bandung 3589 3753120 38
## 2 2 Bandung Barat 3754 1884190 18
## 3 3 Bekasi 1902 3273870 7
## 4 4 Bogor 3404 5682300 23
## 5 5 Ciamis 1420 1259230 11
## 6 6 Cianjur 1932 2584990 11summary(data$JP)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 209790 1064345 1884190 1864637 2569685 5682300data$kab_kota## [1] "Bandung" "Bandung Barat" "Bekasi" "Bogor"
## [5] "Ciamis" "Cianjur" "Cirebon" "Garut"
## [9] "Indramayu" "Karawang" "Kota Bandung" "Kota Banjar"
## [13] "Kota Bekasi" "Kota Bogor" "Kota Cimahi" "Kota Cirebon"
## [17] "Kota Depok" "Kota Sukabumi" "Kota Tasikmalaya" "Kuningan"
## [21] "Majalengka" "Pangandaran" "Purwakarta" "Subang"
## [25] "Sukabumi" "Sumedang" "Tasikmalaya"setdiff(data$kab_kota, Indo$WADMKK)## character(0)# Merged data
jabar_merged <- Indo %>%
left_join(data, by = c("WADMKK" = "kab_kota"))# ------------------------------- HISTOGRAM DAN BOXPLOT --------------------------------
library(openxlsx)
library(ggplot2)
# -------------------------------
# Variabel: JK
# -------------------------------
# Tampilkan ringkasan tanpa variabel no dan kab_kota
summary(data[, !(names(data) %in% c("no", "kab_kota"))])## JK JP KDBD
## Min. : 400 Min. : 209790 Min. : 0.00
## 1st Qu.: 989 1st Qu.:1064345 1st Qu.: 5.50
## Median :1902 Median :1884190 Median :11.00
## Mean :2275 Mean :1864637 Mean :12.56
## 3rd Qu.:3274 3rd Qu.:2569685 3rd Qu.:17.00
## Max. :7680 Max. :5682300 Max. :38.00ggplot(data, aes(x = JK)) +
geom_histogram(fill = "skyblue", color = "black", bins = 15) +
theme_minimal() +
labs(title = "Histogram Jumlah Kasus DBD (JK)", x = "Jumlah Kasus", y = "Frekuensi")
ggplot(data, aes(y = JK)) +
geom_boxplot(fill = "orange", color = "black") +
theme_minimal() +
labs(title = "Boxplot Jumlah Kasus DBD (JK)", y = "Jumlah Kasus")
# ---------------------------- LIBRARY ----------------------------
library(sf)
library(ggplot2)
library(dplyr)
# ---------------------------- KLASIFIKASI JUMLAH KASUS ----------------------------
# Buat kategori jumlah kasus DBD
breaks <- c(0, 1000, 2000, 3000, Inf)
labels <- c("0–1000", "1001–2000", "2001–3000", ">3000")
jabar_merged$Kategori_JK <- cut(jabar_merged$JK,
breaks = breaks,
labels = labels,
right = TRUE)
# ---------------------------- VISUALISASI PETA SEBARAN ----------------------------
ggplot() +
geom_sf(data = jabar_merged, aes(fill = Kategori_JK), color = "white", size = 0.25) +
scale_fill_manual(values = c("0–1000" = "#ffffb2",
"1001–2000" = "#fecc5c",
"2001–3000" = "#fd8d3c",
">3000" = "#e31a1c")) +
labs(title = "Peta Sebaran Jumlah Kasus DBD di Jawa Barat",
fill = "Jumlah Kasus (JK)") +
theme_minimal() +
theme(axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
# -------------------------------
# Variabel: JP
# -------------------------------
# Tampilkan ringkasan tanpa variabel no dan kab_kota
summary(data[, !(names(data) %in% c("no", "kab_kota"))])## JK JP KDBD
## Min. : 400 Min. : 209790 Min. : 0.00
## 1st Qu.: 989 1st Qu.:1064345 1st Qu.: 5.50
## Median :1902 Median :1884190 Median :11.00
## Mean :2275 Mean :1864637 Mean :12.56
## 3rd Qu.:3274 3rd Qu.:2569685 3rd Qu.:17.00
## Max. :7680 Max. :5682300 Max. :38.00ggplot(data, aes(x = JP)) +
geom_histogram(fill = "skyblue", color = "black", bins = 15) +
theme_minimal() +
labs(title = "Histogram Jumlah Penduduk (JP)", x = "Jiwa", y = "Frekuensi")
ggplot(data, aes(y = JP)) +
geom_boxplot(fill = "orange", color = "black") +
theme_minimal() +
labs(title = "Boxplot Jumlah Penduduk (JP)", y = "Jiwa")
# ---------------------------- KLASIFIKASI JUMLAH KASUS ----------------------------
# Buat kategori jumlah kasus DBD
breaks <- c(209000, 1064000, 1884000, 2570000, Inf)
labels <- c("209000–1064000", "1064001–1884000", "1884001–2570000", ">2570000")
jabar_merged$Kategori_JP <- cut(jabar_merged$JP,
breaks = breaks,
labels = labels,
right = TRUE)
# ---------------------------- VISUALISASI PETA SEBARAN ----------------------------
ggplot() +
geom_sf(data = jabar_merged, aes(fill = Kategori_JP), color = "white", size = 0.25) +
scale_fill_manual(values = c("209000–1064000" = "#ffffb2",
"1064001–1884000" = "#fecc5c",
"1884001–2570000" = "#fd8d3c",
">2570000" = "#e31a1c")) +
labs(title = "Peta Sebaran Jumlah Penduduk di Jawa Barat",
fill = "Jumlah Penduduk (JP)") +
theme_minimal() +
theme(axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
# ---------------------------- HITUNG CFR (%) ----------------------------
jabar_merged <- jabar_merged %>%
mutate(CFR = (KDBD / JK) * 100)
# Cek ringkasan untuk menentukan batas
summary(jabar_merged$CFR)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.3515 0.5133 0.6104 0.7440 1.4118# ---------------------------- KLASIFIKASI CFR DIBULATKAN (DENGAN NILAI 0 MASUK KELAS TERENDAH) ----------------------------
breaks_cfr <- c(0, 0.35, 0.51, 0.74, Inf)
labels_cfr <- c("≤0.35%", "0.36–0.51%", "0.52–0.74%", ">0.74%")
jabar_merged$Kategori_CFR <- cut(
jabar_merged$CFR,
breaks = breaks_cfr,
labels = labels_cfr,
right = TRUE,
include.lowest = TRUE
)
# ---------------------------- VISUALISASI PETA CFR (%) ----------------------------
ggplot() +
geom_sf(data = jabar_merged, aes(fill = Kategori_CFR), color = "black", size = 0.25) +
scale_fill_manual(values = c("≤0.35%" = "#ffffb2",
"0.36–0.51%" = "#fecc5c",
"0.52–0.74%" = "#fd8d3c",
">0.74%" = "#e31a1c")) +
labs(title = "Peta Sebaran Case Fatality Rate (CFR) DBD di Jawa Barat (%)",
fill = "CFR (%)") +
theme_minimal() +
theme(axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
# ---------------------------- PETA DISTRIBUSI PENYAKIT: ANALISIS SPASIAL ----------------------------
## ---------------------------- PETA CHOROPLETH DISTRIBUSI DBD ----------------------------
### ---------------------------- HITUNG PREVALENSI (%) ----------------------------
jabar_merged <- jabar_merged %>%
mutate(Prevalensi = (JK / JP) * 100)
# Cek ringkasan untuk menentukan batas kelas
summary(jabar_merged$Prevalensi)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02759 0.08551 0.12837 0.15133 0.19662 0.44595# ---------------------------- KLASIFIKASI PREVALENSI DIBULATKAN ----------------------------
breaks_prev <- c(0, 0.09, 0.13, 0.20, Inf)
labels_prev <- c("≤0.09%", "0.10–0.13%", "0.14–0.20%", ">0.20%")
jabar_merged$Kategori_Prev <- cut(
jabar_merged$Prevalensi,
breaks = breaks_prev,
labels = labels_prev,
right = TRUE,
include.lowest = TRUE
)
# ---------------------------- VISUALISASI PETA PREVALENSI (%) ----------------------------
ggplot() +
geom_sf(data = jabar_merged, aes(fill = Kategori_Prev), color = "black", size = 0.25) +
scale_fill_manual(values = c("≤0.09%" = "#ffffb2",
"0.10–0.13%" = "#fecc5c",
"0.14–0.20%" = "#fd8d3c",
">0.20%" = "#e31a1c")) +
labs(title = "Peta Sebaran Prevalensi DBD di Jawa Barat (%)",
fill = "Prevalensi (%)") +
theme_minimal() +
theme(axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
## ---------------------------- PETA GLOBAL MORAN'S I ----------------------------
# Untuk mengetahui apakah penyebaran DBD membentuk cluster atau acak.
library(spdep)
# Buat neighbor dan weight matrix
neighbors <- poly2nb(jabar_merged)
weights <- nb2listw(neighbors, style = "W")
# Uji Global Moran's I
moran_result <- moran.test(jabar_merged$Prevalensi, weights)
moran_result##
## Moran I test under randomisation
##
## data: jabar_merged$Prevalensi
## weights: weights
##
## Moran I statistic standard deviate = -2.3287, p-value = 0.9901
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.35072320 -0.03846154 0.01798106## ---------------------------- Analisis Cluster Lokal (LISA Map – Local Moran’s I) ----------------------------
local_moran <- localmoran(jabar_merged$Prevalensi, weights)
# Tambahkan hasil ke data
jabar_merged$lisa_I <- local_moran[, 1]
jabar_merged$lisa_p <- local_moran[, 5]
# Klasifikasi cluster
jabar_merged$lisa_cluster <- case_when(
jabar_merged$lisa_p < 0.05 & jabar_merged$lisa_I > 0 ~ "High-High",
jabar_merged$lisa_p < 0.05 & jabar_merged$lisa_I < 0 ~ "Low-Low",
TRUE ~ "Non-Significant"
)
# Peta LISA
ggplot(jabar_merged) +
geom_sf(aes(fill = lisa_cluster), color = "black", size = 0.3) +
scale_fill_manual(values = c(
"High-High" = "red",
"Low-Low" = "blue",
"Non-Significant" = "gray80"
)) +
labs(
title = "Peta LISA Kasus DBD Jawa Barat 2024",
fill = "Tipe Cluster"
) +
theme_minimal()
## ---------------------------- Analisis Hotspot Getis–Ord Gi* ----------------------------
gi_star <- localG(jabar_merged$Prevalensi, weights)
jabar_merged$gi_z <- as.numeric(gi_star)
ggplot(jabar_merged) +
geom_sf(aes(fill = gi_z), color = "gray30", size = 0.3) +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0, name = "Z-score (Gi*)") +
labs(title = "Hotspot Kasus DBD (Getis–Ord Gi*) Jawa Barat 2024") +
theme_minimal()