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 TRUE
library(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, union
library(ggplot2)
library(openxlsx)
## Warning: package 'openxlsx' was built under R version 4.3.3
library(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.3
library(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/Spasial/RBI_50K_2023_Jawa Barat.x26272/RBI_50K_2023_Jawa Barat.shp")
## Reading layer `RBI_50K_2023_Jawa Barat' from data source 
##   `D:\Kuliah\Sem 5\Spasial\RBI_50K_2023_Jawa Barat.x26272\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 84
names(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 --------------------------------
data <- read.xlsx("D:/Kuliah/Sem 5/Epidem/laporan epi/UAS/TBC fix.xlsx")
head(data)
##        kab_kota    JK      JP HIV   ST    DB Sanitasi Faskes Miskin   KP
## 1       Bandung 14012 3753120 598 8.87 44200    96.79     79  239.9 2156
## 2 Bandung Barat  4569 1884190 153 4.95 15911    94.91     38  179.7 1468
## 3        Bekasi 15458 3273870 893 1.70 38666    97.00    102  204.5 2617
## 4         Bogor 29110 5682300 831 1.91 65751    88.77    129  446.8 1899
## 5        Ciamis  3477 1259230 110 2.52 14906    95.16     42   90.8  789
## 6       Cianjur 12197 2584990 293 2.51 15838    69.99     54  239.3  712
##    Rate_1k
## 1 3.733427
## 2 2.424915
## 3 4.721629
## 4 5.122926
## 5 2.761211
## 6 4.718393
summary(data$JK)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     968    4318    6092    8507   12304   29110
data$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 --------------------------------
# -------------------------------
# Variabel: JK
# -------------------------------

# Tampilkan ringkasan tanpa variabel no dan kab_kota
summary(data[, !(names(data) %in% c("no", "kab_kota"))])
##        JK              JP               HIV             ST        
##  Min.   :  968   Min.   : 209790   Min.   :  51   Min.   : 1.600  
##  1st Qu.: 4318   1st Qu.:1064345   1st Qu.: 158   1st Qu.: 2.405  
##  Median : 6092   Median :1884190   Median : 295   Median : 3.580  
##  Mean   : 8507   Mean   :1864637   Mean   : 394   Mean   : 5.437  
##  3rd Qu.:12304   3rd Qu.:2569685   3rd Qu.: 492   3rd Qu.: 8.085  
##  Max.   :29110   Max.   :5682300   Max.   :1400   Max.   :17.080  
##        DB           Sanitasi          Faskes           Miskin     
##  Min.   : 2415   Min.   : 34.56   Min.   : 14.00   Min.   : 11.2  
##  1st Qu.:13425   1st Qu.: 83.67   1st Qu.: 34.00   1st Qu.: 75.3  
##  Median :17912   Median : 95.78   Median : 48.00   Median :131.8  
##  Mean   :24749   Mean   : 89.45   Mean   : 54.74   Mean   :142.5  
##  3rd Qu.:38736   3rd Qu.: 98.55   3rd Qu.: 73.00   3rd Qu.:196.2  
##  Max.   :66896   Max.   :100.00   Max.   :129.00   Max.   :446.8  
##        KP             Rate_1k      
##  Min.   :  385.0   Min.   : 1.229  
##  1st Qu.:  832.5   1st Qu.: 3.321  
##  Median : 1468.0   Median : 4.015  
##  Mean   : 3910.9   Mean   : 5.844  
##  3rd Qu.: 5826.0   3rd Qu.: 5.888  
##  Max.   :15176.0   Max.   :33.934
ggplot(data, aes(x = JK)) +
  geom_histogram(fill = "skyblue", color = "black", bins = 15) +
  theme_minimal() +
  labs(title = "Histogram Jumlah Kasus TBC (JK)", x = "Jumlah Kasus", y = "Frekuensi")

ggplot(data, aes(y = JK)) +
  geom_boxplot(fill = "orange", color = "black") +
  theme_minimal() +
  labs(title = "Boxplot Jumlah Kasus TBC (JK)", y = "Jumlah Kasus")

# ---------------------------- KLASIFIKASI JUMLAH KASUS ----------------------------
# Buat kategori jumlah kasus TBC
breaks <- c(0, 4300, 6100, 12300, Inf)
labels <- c("0–4300", "4301–6100", "6101–12300", ">12300")

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–4300" = "#ffffb2",
                               "4301–6100" = "#fecc5c",
                               "6101–12300" = "#fd8d3c",
                               ">12300" = "#e31a1c")) +
  labs(title = "Peta Sebaran Jumlah Kasus TBC 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               HIV             ST        
##  Min.   :  968   Min.   : 209790   Min.   :  51   Min.   : 1.600  
##  1st Qu.: 4318   1st Qu.:1064345   1st Qu.: 158   1st Qu.: 2.405  
##  Median : 6092   Median :1884190   Median : 295   Median : 3.580  
##  Mean   : 8507   Mean   :1864637   Mean   : 394   Mean   : 5.437  
##  3rd Qu.:12304   3rd Qu.:2569685   3rd Qu.: 492   3rd Qu.: 8.085  
##  Max.   :29110   Max.   :5682300   Max.   :1400   Max.   :17.080  
##        DB           Sanitasi          Faskes           Miskin     
##  Min.   : 2415   Min.   : 34.56   Min.   : 14.00   Min.   : 11.2  
##  1st Qu.:13425   1st Qu.: 83.67   1st Qu.: 34.00   1st Qu.: 75.3  
##  Median :17912   Median : 95.78   Median : 48.00   Median :131.8  
##  Mean   :24749   Mean   : 89.45   Mean   : 54.74   Mean   :142.5  
##  3rd Qu.:38736   3rd Qu.: 98.55   3rd Qu.: 73.00   3rd Qu.:196.2  
##  Max.   :66896   Max.   :100.00   Max.   :129.00   Max.   :446.8  
##        KP             Rate_1k      
##  Min.   :  385.0   Min.   : 1.229  
##  1st Qu.:  832.5   1st Qu.: 3.321  
##  Median : 1468.0   Median : 4.015  
##  Mean   : 3910.9   Mean   : 5.844  
##  3rd Qu.: 5826.0   3rd Qu.: 5.888  
##  Max.   :15176.0   Max.   :33.934
ggplot(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 PENDUDUK ----------------------------

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))

# ----------------------------  Ukuran Epidem Jabar ----------------------------

#Total_JP & Total_JK 
total_JP <- sum(data$JP, na.rm = TRUE)
total_JK <- sum(data$JK, na.rm = TRUE)

total_JP
## [1] 50345190
total_JK
## [1] 229683
# Prevalensi_1k 
prev_JB_1k <- (total_JK / total_JP) * 1000
prev_JB_1k
## [1] 4.562164
# Relative Risk (RR) per kabupaten/kota
data$RR <- (data$JK / data$JP) / (total_JK / total_JP)
head(data[, c("kab_kota", "Rate_1k", "RR")])
##        kab_kota  Rate_1k        RR
## 1       Bandung 3.733427 0.8183457
## 2 Bandung Barat 2.424915 0.5315273
## 3        Bekasi 4.721629 1.0349539
## 4         Bogor 5.122926 1.1229158
## 5        Ciamis 2.761211 0.6052416
## 6       Cianjur 4.718393 1.0342447
summary(data$RR)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2695  0.7280  0.8800  1.2809  1.2907  7.4381
# ---------------------------- HITUNG RR/ Relative Risk ----------------------------

# Rata-rata Resiko Jawa Barat (Total Kasus/Total Penduduk)

risk_jabar <- sum(jabar_merged$JK, na.rm = TRUE) / 
  sum(jabar_merged$JP, na.rm = TRUE)

jabar_merged <- jabar_merged %>%
  mutate(RR = (JK / JP) / risk_jabar)

summary(jabar_merged$RR)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2695  0.7280  0.8800  1.2809  1.2907  7.4381
# ---------------------------- KLASIFIKASI RR DIBULATKAN (DENGAN NILAI 0 MASUK KELAS TERENDAH) ----------------------------
breaks_rr <- c(0, 0.70, 0.90, 1.30, Inf)
labels_rr <- c("≤0.70", "0.71–0.90", "0.91–1.30", ">1.30")

jabar_merged$Kategori_RR <- cut(
  jabar_merged$RR,
  breaks = breaks_rr,
  labels = labels_rr,
  right = TRUE,
  include.lowest = TRUE
)

# ---------------------------- VISUALISASI PETA RR  ----------------------------
ggplot() +
  geom_sf(data = jabar_merged, aes(fill = Kategori_RR), color = "black", size = 0.25) +
  scale_fill_manual(values = c("≤0.70" = "#ffffb2",
                               "0.71–0.90" = "#fecc5c",
                               "0.91–1.30" = "#fd8d3c",
                               ">1.30" = "#e31a1c")) +
  labs(title = "Peta Sebaran Relative Risk (RR) TBC di Jawa Barat",
       fill = "RR") +
  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 TBC ----------------------------

###  ---------------------------- HITUNG PREVALENSI ----------------------------
jabar_merged <- jabar_merged %>%
  mutate(Prevalensi_1k = (JK / JP) * 1000)

# Cek ringkasan untuk menentukan batas kelas
summary(jabar_merged$Prevalensi_1k)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.229   3.321   4.015   5.844   5.888  33.934
# ---------------------------- KLASIFIKASI PREVALENSI DIBULATKAN ----------------------------
breaks_prev <- c(0, 3.3, 4.0, 6.0, Inf)
labels_prev <- c("≤3.3", "3.4–4.0", "4.1–6.0", ">6.0")

jabar_merged$Kategori_Prev <- cut(
  jabar_merged$Prevalensi_1k,
  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("≤3.3" = "#ffffb2",
                               "3.4–4.0" = "#fecc5c",
                               "4.1–6.0" = "#fd8d3c",
                               ">6.0" = "#e31a1c")) +
  labs(title = "Peta Sebaran Prevalensi TBC di Jawa Barat",
       fill = "Prevalensi per 1000 penduduk") +
  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 TBC 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_1k, weights)
moran_result
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$Prevalensi_1k  
## weights: weights    
## 
## Moran I statistic standard deviate = -1.5969, p-value = 0.9449
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.171094749      -0.038461538       0.006898836
## ---------------------------- Analisis Cluster Lokal (LISA Map – Local Moran’s I) ----------------------------

local_moran <- localmoran(jabar_merged$Prevalensi_1k, 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 TBC Jawa Barat 2024",
    fill = "Tipe Cluster"
  ) +
  theme_minimal()

## ---------------------------- Analisis Hotspot Getis–Ord Gi* ----------------------------

gi_star <- localG(jabar_merged$Prevalensi_1k, 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 TBC (Getis–Ord Gi*) Jawa Barat 2024") +
  theme_minimal()

# ---------------------------- PEMODELAN PEMETAAN ----------------------------
# Paket yang dibutuhkan
library(car)        # VIF
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(lmtest)     # bptest
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(tseries)    # jarque.bera.test
## Warning: package 'tseries' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(MASS)       # boxcox (opsional)
## Warning: package 'MASS' was built under R version 4.3.3
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(ggplot2)

# Pastikan data berbentuk data.frame
str(jabar_merged)
## Classes 'sf' and 'data.frame':   27 obs. of  46 variables:
##  $ NAMOBJ       : chr  "Bandung" "Bandung Barat" "Bekasi" "Bogor" ...
##  $ FCODE        : chr  "BA03050040" "BA03050040" "BA03050040" "BA03050040" ...
##  $ REMARK       : chr  NA NA NA NA ...
##  $ METADATA     : chr  "TASWIL5000020230907KABKOTA" "TASWIL5000020230907KABKOTA" "TASWIL5000020230907KABKOTA" "TASWIL5000020230907KABKOTA" ...
##  $ SRS_ID       : chr  "4326" "4326" "4326" "4326" ...
##  $ KDBBPS       : chr  NA NA NA NA ...
##  $ KDCBPS       : chr  NA NA NA NA ...
##  $ KDCPUM       : chr  NA NA NA NA ...
##  $ KDEBPS       : chr  NA NA NA NA ...
##  $ KDEPUM       : chr  NA NA NA NA ...
##  $ KDPBPS       : chr  NA NA NA NA ...
##  $ KDPKAB       : chr  "32.04" "32.17" "32.16" "32.01" ...
##  $ KDPPUM       : chr  "32" "32" "32" "32" ...
##  $ LUASWH       : num  1741 1283 1251 2992 1596 ...
##  $ TIPADM       : num  4 4 4 4 4 4 4 4 4 4 ...
##  $ WADMKC       : chr  NA NA NA NA ...
##  $ WADMKD       : chr  NA NA NA NA ...
##  $ WADMKK       : chr  "Bandung" "Bandung Barat" "Bekasi" "Bogor" ...
##  $ WADMPR       : chr  "Jawa Barat" "Jawa Barat" "Jawa Barat" "Jawa Barat" ...
##  $ WIADKC       : chr  NA NA NA NA ...
##  $ WIADKK       : chr  NA "Bandung" NA NA ...
##  $ WIADPR       : chr  NA NA NA NA ...
##  $ WIADKD       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ SHAPE_Leng   : num  3.24 2.56 3.14 4.53 3 ...
##  $ SHAPE_Area   : num  0.142 0.105 0.102 0.245 0.131 ...
##  $ JK           : num  14012 4569 15458 29110 3477 ...
##  $ JP           : num  3753120 1884190 3273870 5682300 1259230 ...
##  $ HIV          : num  598 153 893 831 110 293 400 328 573 886 ...
##  $ ST           : num  8.87 4.95 1.7 1.91 2.52 ...
##  $ DB           : num  44200 15911 38666 65751 14906 ...
##  $ Sanitasi     : num  96.8 94.9 97 88.8 95.2 ...
##  $ Faskes       : num  79 38 102 129 42 54 71 76 60 75 ...
##  $ Miskin       : num  239.9 179.7 204.5 446.8 90.8 ...
##  $ KP           : num  2156 1468 2617 1899 789 ...
##  $ Rate_1k      : num  3.73 2.42 4.72 5.12 2.76 ...
##  $ geometry     :sfc_MULTIPOLYGON of length 27; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:13089, 1:2] 108 108 108 108 108 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  $ Kategori_JK  : Factor w/ 4 levels "0–4300","4301–6100",..: 4 2 4 4 1 3 3 3 2 4 ...
##  $ Kategori_JP  : Factor w/ 4 levels "209000–1064000",..: 4 3 4 4 2 4 3 4 3 3 ...
##  $ RR           : num  0.818 0.532 1.035 1.123 0.605 ...
##  $ Kategori_RR  : Factor w/ 4 levels "≤0.70","0.71–0.90",..: 2 1 3 3 1 3 2 2 1 3 ...
##  $ Prevalensi_1k: num  3.73 2.42 4.72 5.12 2.76 ...
##  $ Kategori_Prev: Factor w/ 4 levels "≤3.3","3.4–4.0",..: 2 1 3 3 1 3 3 2 1 3 ...
##  $ lisa_I       : num  0.0687 0.0542 0.0204 -0.0648 0.1417 ...
##  $ lisa_p       : num  0.5321 0.7512 0.8362 0.0796 0.4238 ...
##  $ lisa_cluster : chr  "Non-Significant" "Non-Significant" "Non-Significant" "Non-Significant" ...
##  $ gi_z         : num  -0.625 -0.317 -0.207 1.753 -0.8 ...
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:45] "NAMOBJ" "FCODE" "REMARK" "METADATA" ...
# membentuk Model Regresi Linear Berganda 
model_lm <- lm(
  Prevalensi_1k ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
  data = jabar_merged
)

summary(model_lm)
## 
## Call:
## lm(formula = Prevalensi_1k ~ HIV + ST + DB + Sanitasi + Faskes + 
##     Miskin + KP, data = jabar_merged)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1652 -2.5544 -1.7108  0.9932 23.7152 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.4049230 12.0964367  -0.612    0.548
## HIV          0.0023243  0.0112125   0.207    0.838
## ST           0.3601374  0.3917980   0.919    0.370
## DB          -0.0001919  0.0001470  -1.305    0.207
## Sanitasi     0.1044568  0.1146017   0.911    0.373
## Faskes       0.0635273  0.1581795   0.402    0.692
## Miskin       0.0104651  0.0364971   0.287    0.777
## KP           0.0002075  0.0004492   0.462    0.649
## 
## Residual standard error: 6.544 on 19 degrees of freedom
## Multiple R-squared:  0.167,  Adjusted R-squared:  -0.14 
## F-statistic: 0.544 on 7 and 19 DF,  p-value: 0.7907
# Uji Asumsi 
#Linearitas
plot(model_lm, which = 1)

#ramsey reset test
library(lmtest)

resettest(model_lm, power = 2:3, type = "fitted")
## 
##  RESET test
## 
## data:  model_lm
## RESET = 4.1642, df1 = 2, df2 = 17, p-value = 0.03374
# uji normalitas residual 
jarque.bera.test(residuals(model_lm))
## 
##  Jarque Bera Test
## 
## data:  residuals(model_lm)
## X-squared = 153.97, df = 2, p-value < 2.2e-16
# Uji Homoskedastisitas
bptest(model_lm)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_lm
## BP = 4.863, df = 7, p-value = 0.6767
# Uji Multikolinearitas
vif(model_lm)
##       HIV        ST        DB  Sanitasi    Faskes    Miskin        KP 
##  8.088758  1.524629  4.010534  1.769091 13.208417  7.618261  2.670133
###### EVALUASI MODEL 

# Koefisien Determinasi
summary(model_lm)$r.squared
## [1] 0.1669584
summary(model_lm)$adj.r.squared
## [1] -0.1399517
# RMSE
rmse <- sqrt(mean(residuals(model_lm)^2))
rmse
## [1] 5.489385
# Uji F (Simultan )
summary(model_lm)
## 
## Call:
## lm(formula = Prevalensi_1k ~ HIV + ST + DB + Sanitasi + Faskes + 
##     Miskin + KP, data = jabar_merged)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1652 -2.5544 -1.7108  0.9932 23.7152 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.4049230 12.0964367  -0.612    0.548
## HIV          0.0023243  0.0112125   0.207    0.838
## ST           0.3601374  0.3917980   0.919    0.370
## DB          -0.0001919  0.0001470  -1.305    0.207
## Sanitasi     0.1044568  0.1146017   0.911    0.373
## Faskes       0.0635273  0.1581795   0.402    0.692
## Miskin       0.0104651  0.0364971   0.287    0.777
## KP           0.0002075  0.0004492   0.462    0.649
## 
## Residual standard error: 6.544 on 19 degrees of freedom
## Multiple R-squared:  0.167,  Adjusted R-squared:  -0.14 
## F-statistic: 0.544 on 7 and 19 DF,  p-value: 0.7907
######## Metode POisson #########

# 1. Membentuk model regresi Poisson
# Model Regresi Poisson
model_pois <- glm(
  JK ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP +
    offset(log(JP)),
  family = poisson(link = "log"),
  data = jabar_merged
)

summary(model_pois)
## 
## Call:
## glm(formula = JK ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + 
##     KP + offset(log(JP)), family = poisson(link = "log"), data = jabar_merged)
## 
## Coefficients:
##               Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -5.556e+00  2.146e-02 -258.884  < 2e-16 ***
## HIV          1.885e-04  1.791e-05   10.522  < 2e-16 ***
## ST          -1.046e-02  6.880e-04  -15.209  < 2e-16 ***
## DB          -1.340e-05  2.566e-07  -52.241  < 2e-16 ***
## Sanitasi    -1.736e-03  2.043e-04   -8.495  < 2e-16 ***
## Faskes       7.409e-03  2.773e-04   26.715  < 2e-16 ***
## Miskin       4.458e-04  6.160e-05    7.237 4.59e-13 ***
## KP           2.289e-05  8.234e-07   27.795  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 63934  on 26  degrees of freedom
## Residual deviance: 53031  on 19  degrees of freedom
## AIC: 53334
## 
## Number of Fisher Scoring iterations: 5
# 2. Interpretasi koefisien (IRR)
# Incidence Rate Ratio (IRR)
exp(coef(model_pois))
## (Intercept)         HIV          ST          DB    Sanitasi      Faskes 
## 0.003862953 1.000188484 0.989590752 0.999986596 0.998265991 1.007436722 
##      Miskin          KP 
## 1.000445898 1.000022888
# Interval kepercayaan IRR
exp(confint(model_pois))
## Waiting for profiling to be done...
##                  2.5 %      97.5 %
## (Intercept) 0.00370372 0.004028805
## HIV         1.00015338 1.000223605
## ST          0.98825687 0.990925778
## DB          0.99998609 0.999987098
## Sanitasi    0.99786660 0.998666026
## Faskes      1.00688937 1.007984611
## Miskin      1.00032515 1.000566731
## KP          1.00002127 1.000024502
# 3. Uji Overdispersion
# 3.1 Metode Rasio Deviance
# Rasio deviance
deviance(model_pois) / df.residual(model_pois)
## [1] 2791.089
# 3.2 Metode Pearson Chi-Square 
# Rasio Pearson Chi-square
sum(residuals(model_pois, type = "pearson")^2) /
  df.residual(model_pois)
## [1] 4995.069
# 3.3. Uji Formal Overdispersion
library(AER)
## Loading required package: sandwich
## Warning: package 'sandwich' was built under R version 4.3.3
## Loading required package: survival
## Warning: package 'survival' was built under R version 4.3.3
dispersiontest(model_pois)
## 
##  Overdispersion test
## 
## data:  model_pois
## z = 1.3878, p-value = 0.0826
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   3514.702
# 4. EVALUASI MODEL POISSON
# AIC
AIC(model_pois)
## [1] 53333.73
# Log-likelihood
logLik(model_pois)
## 'log Lik.' -26658.87 (df=8)
# 5. DIAGNOSTIK MODEL POISSON
# Plot residual
plot(model_pois)

######## Metode POisson-Gamma #########
#1. Model
library(MASS)
model_pg <- glm.nb(
  JK ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP +
    offset(log(JP)),
  data = jabar_merged
)

# Ringkasan Model
summary(model_pg)
## 
## Call:
## glm.nb(formula = JK ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + 
##     KP + offset(log(JP)), data = jabar_merged, init.theta = 3.029303181, 
##     link = log)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -7.107e+00  1.062e+00  -6.690 2.23e-11 ***
## HIV          6.229e-04  9.846e-04   0.633   0.5270    
## ST           4.113e-02  3.441e-02   1.195   0.2320    
## DB          -3.067e-05  1.291e-05  -2.375   0.0175 *  
## Sanitasi     1.446e-02  1.006e-02   1.437   0.1507    
## Faskes       4.619e-03  1.389e-02   0.332   0.7395    
## Miskin       2.983e-03  3.205e-03   0.931   0.3520    
## KP           5.371e-05  3.946e-05   1.361   0.1734    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.0293) family taken to be 1)
## 
##     Null deviance: 40.390  on 26  degrees of freedom
## Residual deviance: 28.468  on 19  degrees of freedom
## AIC: 542.5
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.029 
##           Std. Err.:  0.784 
## 
##  2 x log-likelihood:  -524.498
#IRR
exp(coef(model_pg))
##  (Intercept)          HIV           ST           DB     Sanitasi       Faskes 
## 0.0008190819 1.0006230738 1.0419824320 0.9999693296 1.0145682216 1.0046296352 
##       Miskin           KP 
## 1.0029878621 1.0000537094
#CI IRR
exp(confint(model_pg))
## Waiting for profiling to be done...
##                    2.5 %      97.5 %
## (Intercept) 9.725871e-05 0.008214414
## HIV         9.984804e-01 1.002709063
## ST          9.754882e-01 1.116204121
## DB          9.999448e-01 0.999995695
## Sanitasi    9.917902e-01 1.035657459
## Faskes      9.764162e-01 1.034328682
## Miskin      9.956450e-01 1.010119624
## KP          9.999744e-01 1.000135530
# AIC & Log-likelihood
AIC(model_pois, model_pg)
##            df        AIC
## model_pois  8 53333.7325
## model_pg    9   542.4983
logLik(model_pois)
## 'log Lik.' -26658.87 (df=8)
logLik(model_pg)
## 'log Lik.' -262.2491 (df=9)
deviance(model_pg) / df.residual(model_pg)
## [1] 1.49833
sum(residuals(model_pg, type = "pearson")^2) /
  df.residual(model_pg)
## [1] 2.003193
#memnaggil INLA

library(INLA)
## Loading required package: Matrix
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.3.3
## This is INLA_23.09.09 built 2023-10-16 17:29:11 UTC.
##  - See www.r-inla.org/contact-us for how to get help.
#################### BYM ########## 

jabar_merged$region_id <- as.integer(factor(jabar_merged$WADMKK))
summary(jabar_merged$region_id)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     7.5    14.0    14.0    20.5    27.0
length(unique(jabar_merged$region_id))
## [1] 27
library(spdep)
library(INLA)

# Pastikan urutan sama
jabar_merged <- jabar_merged[order(jabar_merged$region_id), ]

# Buat neighbors
nb <- poly2nb(jabar_merged)

# Simpan adjacency
nb2INLA("jabar.adj", nb)
file.exists("jabar.adj")
## [1] TRUE
g <- inla.read.graph("jabar.adj")
g$n
## [1] 27
# expected cases (offset)
rate_global <- sum(jabar_merged$JK) / sum(jabar_merged$JP)
jabar_merged$E <- jabar_merged$JP * rate_global

# jalankan model BYM
formula_bym <- JK ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP +
  f(region_id, model = "bym", graph = "jabar.adj") +
  offset(log(E))

model_bym <- inla(
  formula_bym,
  family = "poisson",
  data = jabar_merged,
  control.compute = list(dic = TRUE, waic = TRUE),
  control.predictor = list(compute = TRUE)
)

#Output yang diinterpretassi

#ringakasan model
summary(model_bym)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " keep = keep, working.directory = working.directory, 
##    silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = 
##    debug, .parent.frame = .parent.frame)" ) 
## Time used:
##     Pre = 1.36, Running = 0.79, Post = 0.393, Total = 2.55 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.066 1.232     -3.500   -1.066      1.367 -1.066   0
## HIV          0.000 0.001     -0.002    0.000      0.003  0.000   0
## ST           0.021 0.040     -0.058    0.021      0.100  0.021   0
## DB           0.000 0.000      0.000    0.000      0.000  0.000   0
## Sanitasi     0.007 0.012     -0.016    0.007      0.030  0.007   0
## Faskes       0.005 0.016     -0.027    0.005      0.037  0.005   0
## Miskin       0.002 0.004     -0.006    0.002      0.009  0.002   0
## KP           0.000 0.000      0.000    0.000      0.000  0.000   0
## 
## Random effects:
##   Name     Model
##     region_id BYM model
## 
## Model hyperparameters:
##                                                mean       sd 0.025quant
## Precision for region_id (iid component)        2.46    0.767       1.25
## Precision for region_id (spatial component) 2203.59 2421.159     145.87
##                                             0.5quant 0.975quant   mode
## Precision for region_id (iid component)         2.36       4.24   2.18
## Precision for region_id (spatial component)  1443.93    8631.71 397.17
## 
## Deviance Information Criterion (DIC) ...............: 341.03
## Deviance Information Criterion (DIC, saturated) ....: 53.99
## Effective number of parameters .....................: 26.99
## 
## Watanabe-Akaike information criterion (WAIC) ...: 332.77
## Effective number of parameters .................: 13.51
## 
## Marginal log-Likelihood:  -330.49 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
#Resiko Relatif
jabar_merged$RR <- model_bym$summary.fitted.values$mean / jabar_merged$E

#Peta Resiko Relatif 

library(tmap)

tm_shape(jabar_merged) +
  tm_fill("RR", style = "quantile", title = "Risiko Relatif TBC") +
  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' 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>)'

# Evaluasi MOdel Bayesian
model_bym$dic$dic
## [1] 341.0339
model_bym$waic$waic
## [1] 332.7667