Nama: Fatih Zahrani
NPM:
140610230014
Dosen Pengampu: I Gede Nyoman Mindra
Jaya, M.Si., Ph.D.
Mata Kuliah: Statistik Spasial
S1 Statistika FMIPA UNPAD
Petunjuk
Download folder Dashboard dan jalankan syntax app.R dibawah ini
Link file:
https://drive.google.com/drive/folders/1Ngg_Tgcaw2tC-1soexGVPaJbzkuX_col?usp=sharing
Syntax Dashboard
# ============================================================
# app.R ??? Dashboard Analisis TBC Jawa Barat 2024
# Version: 1.0
# ============================================================
# ----------------------------
# Load libraries
# ----------------------------
library(shiny)
library(shinydashboard)
library(shinyWidgets)
library(plotly)
library(leaflet)
library(DT)
library(spdep)
library(sf)
library(sp)
library(openxlsx)
library(tidyr)
library(ggplot2)
library(DiagrammeR)
library(viridis)
library(scales)
library(MASS)
library(dplyr)
library(terra)
library(spatialreg)
library(performance)
library(conflicted)
conflict_prefer("filter", "dplyr") # Contoh jika ada konflik dplyr
conflicts_prefer(shinydashboard::box)
conflicts_prefer(dplyr::select)
conflicts_prefer(plotly::layout)
# ============================================================
# 1. GLOBAL SECTION (Data & Model)
# ============================================================
# --- Load shapefile dan data TBC ---
Indo <- st_read("shapefile") |>
st_zm(drop = TRUE, what = "ZM")
data <- read.xlsx("data/DATA TBC.xlsx")
# --- Merge shapefile dan data ---
jabar_merged <- Indo %>%
left_join(data, by = c("WADMKK" = "kab_kota"))
# ============================================================
# HITUNG PREVALENSI & RELATIVE RISK (RR)
# ============================================================
# Pastikan numerik
jabar_merged <- jabar_merged %>%
mutate(
JK = as.numeric(JK),
JP = as.numeric(JP)
)
# Hitung risk rata-rata Jawa Barat
risk_jabar <- sum(jabar_merged$JK, na.rm = TRUE) /
sum(jabar_merged$JP, na.rm = TRUE)
# Gabungkan & hitung RR
jabar_merged <- jabar_merged %>%
mutate(
Prev_1k = (JK / JP) * 1000,
RR = (JK / JP) / risk_jabar,
Rate_1k = Prev_1k
)
# --- Total untuk Jawa Barat (indikator agregat) ---
total_JK <- sum(jabar_merged$JK, na.rm = TRUE)
total_JP <- sum(jabar_merged$JP, na.rm = TRUE)
prev_jabar_1k <- (total_JK / total_JP) * 1000
RR_jabar <- 1 # referensi (rata-rata provinsi)
# --- Konversi ke Spatial untuk analisis spasial ---
jabar_merged <- st_zm(jabar_merged, drop = TRUE, what = "ZM")
jabar_sp <- as_Spatial(jabar_merged)
row.names(jabar_sp) <- jabar_sp$WADMKK
# --- Spatial weight matrix ---
W <- poly2nb(jabar_sp, row.names = row.names(jabar_sp), queen = TRUE)
WL <- nb2listw(W, style = "W", zero.policy = TRUE)
# Titik centroid kab/kota
centroid_sf <- st_centroid(jabar_merged)
centroid_sp <- as(centroid_sf, "Spatial")
coords <- st_coordinates(centroid_sf)
jabar_merged$x <- coords[,1]
jabar_merged$y <- coords[,2]
# ===============================
# GRID PREDIKSI SPASIAL
# ===============================
bbox <- st_bbox(jabar_merged)
grd <- expand.grid(
x = seq(bbox["xmin"], bbox["xmax"], length.out = 100),
y = seq(bbox["ymin"], bbox["ymax"], length.out = 100)
)
model_trend <- lm(JK ~ x + y, data = jabar_merged)
grd$pred_JK <- predict(model_trend, newdata = as.data.frame(grd))
coordinates(grd) <- ~ x + y
centroid_sp <- as_Spatial(st_centroid(jabar_merged))
proj4string(grd) <- proj4string(centroid_sp)
proj4string(grd) <- proj4string(centroid_sp)
gridded(grd) <- TRUE
# Data untuk interpolasi
df_interp <- data.frame(
x = coordinates(centroid_sp)[,1],
y = coordinates(centroid_sp)[,2],
Rate_1k = jabar_merged$Rate_1k
)
df_trend <- df_interp
coordinates(df_interp) <- ~ x + y
proj4string(df_interp) <- proj4string(centroid_sp)
"x" %in% names(jabar_merged) # TRUE
"y" %in% names(jabar_merged) # TRUE
# IDW
idw_model <- gstat::gstat(
formula = Rate_1k ~ 1,
data = df_interp,
set = list(idp = 2)
)
# PLOT IDW
idw_pred <- predict(idw_model, grd)
idw_raster <- rast(data.frame(coordinates(idw_pred), idw_pred$var1.pred))
library(gstat)
idw_res <- idw(Rate_1k ~ 1, df_interp, newdata = grd)
spplot(idw_res["var1.pred"])
r <- rast(idw_res)
# Trend Surface
grd_df <- as.data.frame(grd)
colnames(grd_df) <- c("x", "y")
trend_pred <- predict(
lm(Rate_1k ~ poly(x, 2) + poly(y, 2), data = df_trend),
newdata = grd_df
)
trend_raster <- rast(grd)
values(trend_raster) <- trend_pred
plot(trend_raster)
#Kriging
vgm_fit <- gstat::fit.variogram(
variogram(Rate_1k ~ 1, df_interp),
vgm("Sph")
)
kriging_pred <- krige(
formula = Rate_1k ~ 1,
locations = df_interp,
newdata = grd,
model = vgm_fit
)
kriging_raster <- rast(data.frame(coordinates(kriging_pred), kriging_pred$var1.pred))
cat("✅ Semua data TBC & RR berhasil dimuat!\n")
# ============================================================
# 2. UI SECTION
# ============================================================
ui <- dashboardPage(
skin = "blue",
dashboardHeader(title = "Dashboard TBC Jawa Barat 2024", titleWidth = 350),
dashboardSidebar(
width = 270,
sidebarMenu(
menuItem("About Dashboard", tabName = "about", icon = icon("info-circle")),
menuItem("Data Overview", tabName = "data_overview", icon = icon("table")),
menuItem("Determinan TBC", tabName = "determinant", icon = icon("project-diagram")),
menuItem("Statistika Deskriptif", tabName = "descriptive", icon = icon("chart-bar")),
menuItem("Sebaran TBC", tabName = "spatial", icon = icon("map-marked-alt")),
menuItem("Ukuran Epidemiologi", tabName = "epidemiology", icon = icon("calculator")),
menuItem("Pemodelan", tabName = "modeling", icon = icon("microscope")),
menuItem("Analisis Spasial", tabName = "spatial_analysis", icon = icon("globe"))
)
),
dashboardBody(
tags$head(
tags$link(rel = "stylesheet", type = "text/css", href = "custom.css")
),
tabItems(
# =====================================================
# 1. ABOUT DASHBOARD
# =====================================================
tabItem(tabName = "about",
fluidRow(
box(width = 12, title = "Dashboard Tuberkulosis (TBC) Jawa Barat 2024",
status = "primary", solidHeader = TRUE,
h4("Tentang Dashboard"),
p("Dashboard ini menyajikan analisis epidemiologi Tuberkulosis (TBC) di Provinsi Jawa Barat tahun 2024.
Fokus utama dashboard adalah analisis sebaran spasial, prevalensi, dan risiko relatif (Relative Risk)
TBC antar kabupaten/kota."),
hr(),
h4("Urgensi"),
tags$ul(
tags$li("TBC masih menjadi salah satu penyakit menular prioritas dengan beban tinggi di Jawa Barat."),
tags$li("Analisis spasial diperlukan untuk mengidentifikasi wilayah dengan risiko relatif lebih tinggi
dibandingkan rata-rata provinsi.")
),
hr(),
h4("Sumber Data"),
p(strong("Kasus TBC, HIV, Diabetes, Sanitasi, Faskes:"), " Open Data Jawa Barat 2024"),
p(strong("Penduduk, Kepadatan, Kemiskinan:"), " BPS Jawa Barat 2024"),
p(strong("Shapefile:"), " RBI 50K Jawa Barat"),
p(strong("Tahun:"), " 2024"),
p(strong("Jumlah Kabupaten/Kota:"), nrow(data)),
hr(),
h4("Penyusun"),
tags$ul(
tags$li("Fatih Zahrani (140610230014)")
),
p(em("Dashboard Version 2.0 - 2025"))
)
)
)
),
# =====================================================
# 2. DATA OVERVIEW
# =====================================================
tabItem(tabName = "data_overview",
fluidRow(
valueBoxOutput("total_cases", width = 4), # Total kasus TBC
valueBoxOutput("max_cases", width = 4),
valueBoxOutput("num_districts", width = 4) # Jumlah wilayah
),
fluidRow(
box(width = 12, title = "Tabel Data TBC", status = "info", solidHeader = TRUE,
sliderInput("n_rows", "Jumlah Baris:", min = 5, max = 50, value = 10),
downloadButton("download_data", "Download Data (Excel)", class = "btn-primary"),
hr(),
DTOutput("data_table"))
)
),
# =====================================================
# 3. DETERMINAN TBC
# =====================================================
tabItem(tabName = "determinant",
fluidRow(
box(width = 12, title = "Faktor Determinan TBC", status = "primary", solidHeader = TRUE,
p("Faktor-faktor sosial, lingkungan, dan kesehatan yang berpotensi memengaruhi kejadian TBC."),
tags$ul(
tags$li("Kasus HIV"),
tags$li("Diabetes Melitus"),
tags$li("Stunting"),
tags$li("Kepadatan Penduduk"),
tags$li("Kemiskinan"),
tags$li("Fasilitas Kesehatan"),
tags$li("Sanitasi Layak")
)
)
)
),
# =====================================================
# 4. STATISTIKA DESKRIPTIF
# =====================================================
tabItem(tabName = "descriptive",
fluidRow(
box(
width = 12,
title = "Statistik Deskriptif Data TBC",
status = "primary",
solidHeader = TRUE,
collapsible = TRUE,
pickerInput(
inputId = "desc_vars",
label = "Pilih Variabel:",
choices = NULL, # ← jangan ditulis manual
multiple = TRUE,
options = list(actionsBox = TRUE, liveSearch = TRUE, selectedTextFormat = "count > 3")
),
tabsetPanel(
tabPanel(
"Histogram",
br(),
plotlyOutput("desc_histogram", height = "500px")
),
tabPanel(
"Boxplot",
br(),
plotlyOutput("desc_boxplot", height = "500px")
),
tabPanel(
"Ringkasan Statistik",
br(),
DTOutput("desc_summary")
)
),
br(),
downloadButton("download_desc", "Download Statistik Deskriptif")
)
)
),
# =====================================================
# 5. SEBARAN TBC
# =====================================================
tabItem(tabName = "spatial",
fluidRow(
valueBoxOutput("total_kasus", width = 4),
valueBoxOutput("max_rr", width = 4),
valueBoxOutput("total_wilayah", width = 4)
),
fluidRow(
box(width = 12, title = "Peta Sebaran Jumlah Kasus TBC",
status = "danger", solidHeader = TRUE,
leafletOutput("map_jk", height = "500px"))
),
fluidRow(
box(
width = 12,
title = "Analisis Autokorelasi Spasial Relative Risk (RR) TBC",
status = "primary",
solidHeader = TRUE,
tabsetPanel(
tabPanel(
"Moran's I & Geary's C",
br(),
fluidRow(
valueBoxOutput("moranI_box", width = 6),
valueBoxOutput("gearyC_box", width = 6)
),
plotlyOutput("moran_plot", height = "500px")
),
tabPanel(
"LISA Cluster Map",
br(),
p("Local Indicators of Spatial Association (LISA) menunjukkan cluster RR TBC."),
leafletOutput("lisa_map", height = "600px")
),
tabPanel(
"Getis-Ord Gi*",
br(),
p("Identifikasi hotspot dan coldspot Relative Risk (RR) TBC."),
leafletOutput("getisord_map", height = "600px")
)
)
)
)
),
# =====================================================
# 6. UKURAN EPIDEMIOLOGI
# =====================================================
tabItem(
tabName = "epidemiology",
# =========================
# VALUE BOX
# =========================
fluidRow(
valueBoxOutput("prev_jabar_box", width = 6),
valueBoxOutput("rr_jabar_box", width = 6)
),
# =========================
# PREVALENSI TBC
# =========================
fluidRow(
box(
width = 12,
title = "Analisis Prevalensi TBC (per 1.000 Penduduk)",
status = "info",
solidHeader = TRUE,
collapsible = TRUE,
tabsetPanel(
tabPanel(
"Peta Prevalensi",
br(),
leafletOutput("map_prevalensi", height = "500px")
),
tabPanel(
"Prevalensi per Wilayah",
br(),
pickerInput(
"prev_cities",
"Pilih Kabupaten/Kota:",
choices = NULL,
multiple = TRUE,
options = list(actionsBox = TRUE, liveSearch = TRUE)
),
fluidRow(
column(6, DTOutput("prev_table")),
column(6, plotlyOutput("prev_bar", height = "400px"))
)
)
)
)
),
# =========================
# RELATIVE RISK (RR)
# =========================
fluidRow(
box(
width = 12,
title = "Analisis Relative Risk (RR) TBC",
status = "warning",
solidHeader = TRUE,
collapsible = TRUE,
tabsetPanel(
tabPanel(
"Peta Relative Risk",
br(),
leafletOutput("map_rr", height = "500px")
),
tabPanel(
"Relative Risk per Wilayah",
br(),
pickerInput(
"rr_cities",
"Pilih Kabupaten/Kota:",
choices = NULL,
multiple = TRUE,
options = list(actionsBox = TRUE, liveSearch = TRUE)
),
fluidRow(
column(6, DTOutput("rr_table")),
column(6, plotlyOutput("rr_bar", height = "400px"))
)
)
)
)
)
),
# =====================================================
# 7. Pemodelan
# =====================================================
tabItem(tabName = "modeling",
fluidRow(
tabBox(
title = "Analisis Pemodelan TBC",
id = "tab_modeling", width = 12,
# --- TAB 1: OLS ---
tabPanel("Regresi OLS (Prevalensi)",
fluidRow(
valueBoxOutput("ols_r2", width = 4),
valueBoxOutput("ols_adj_r2", width = 4),
valueBoxOutput("ols_rmse", width = 4)
),
fluidRow(
column(7, box(title = "Ringkasan Model OLS", width = NULL, status = "primary", verbatimTextOutput("ols_summary"))),
column(5, box(title = "Uji Asumsi Klasik", width = NULL, status = "warning", DTOutput("ols_asumsi_table")))
),
plotOutput("ols_plot")
),
# --- TAB 2: POISSON & NEGATIVE BINOMIAL ---
tabPanel("Poisson-Gamma",
fluidRow(
valueBoxOutput("dispersion_status", width = 6),
valueBoxOutput("best_model_aic", width = 6)
),
hr(),
h4("Tabel 6. Interpretasi Nilai Incident Rate Ratio (IRR) - Poisson-Gamma"),
DTOutput("irr_pg_table"),
hr()
),
# --- TAB 3: BYM-INLA (SPASIAL) ---
tabPanel("Bayesian BYM-INLA",
# 1. Baris ValueBox (DIC & WAIC)
# width = 6 + 6 = 12 (Memenuhi satu baris penuh)
fluidRow(
valueBoxOutput("bym_dic", width = 6),
valueBoxOutput("bym_waic", width = 6)
),
# 2. Baris Peta (FULL WIDTH)
fluidRow(
box(
title = "Peta Risiko Relatif (Posterior RR)",
width = 12, # <--- UBAH JADI 12 (Full Width)
status = "primary", # Opsional: Memberi warna header biru
solidHeader = TRUE, # Opsional: Header berwarna solid
leafletOutput("map_bym_rr", height = "600px") # Tinggi diperbesar agar proporsional
)
),
# 3. Baris Tabel Fixed Effects
# (Karena peta sudah full, tabel ini sebaiknya ditaruh di bawahnya selebar 12 juga)
fluidRow(
box(
title = "Fixed Effects (Estimasi Bayesian)",
width = 12,
status = "info",
DTOutput("bym_fixed_table")
)
),
hr(),
# 4. Tabel Interpretasi CI
h4("Tabel 7. Interval Kepercayaan (CI) IRR dan Interpretasi Konsistensi - BYM"),
DTOutput("bym_ci_table")
)
)
)
),
# =====================================================
# 8. Analisis Spasial
# =====================================================
tabItem(tabName = "spatial_analysis",
fluidRow(
tabBox(
title = "Analisis Spasial Klasik dan Interpolasi",
id = "spatial_analysis", width = 12,
tabsetPanel(
type = "tabs",
tabPanel(
"Model OLS",
h4("Bentuk Matematis Model"),
withMathJax(uiOutput("ols_formula")),
br(),
h4("Ringkasan Model OLS"),
tableOutput("ols_table"),
br(),
h4("Uji Autokorelasi Spasial Residual"),
tableOutput("moran_ols"),
br(),
h4("Interpretasi"),
verbatimTextOutput("ols_interpretation")
),
tabPanel(
"Model Spasial Klasik",
h4("Ringkasan Model Spasial"),
verbatimTextOutput("spatial_classic_table"),
br(),
h4("Uji Dependensi Spasial"),
DTOutput("lm_test"),
br(),
h4("Interpretasi"),
verbatimTextOutput("spatial_classic_interp")
),
tabPanel(
"Interpolasi Spasial",
tabsetPanel(
tabPanel("IDW", plotOutput("idw_plot", height = 500)),
tabPanel("Trend Surface", plotOutput("trend_plot", height = 500)),
tabPanel("Ordinary Kriging", plotOutput("kriging_plot", height = 500))
)
)
)
)
)
)
)
)
# ============================================================
# 3. SERVER SECTION
# ============================================================
server <- function(input, output, session) {
options(error = recover)
# ============================================================
# REACTIVE VALUES
# ============================================================
# Update choices untuk picker inputs
observe({
city_choices <- data$kab_kota
updatePickerInput(session, "desc_cities", choices = city_choices)
updatePickerInput(session, "prev_cities", choices = city_choices)
updatePickerInput(session, "rr_cities", choices = city_choices)
})
# ============================================================
# UPDATE PICKER VARIABEL NUMERIK (STATISTIK DESKRIPTIF)
# ============================================================
observe({
numeric_vars <- names(data)[sapply(data, is.numeric)]
updatePickerInput(
session,
"desc_vars",
choices = numeric_vars,
selected = "JK"
)
})
# ============================================================
# 2. DATA OVERVIEW
# ============================================================
# ValueBoxes
output$total_cases <- renderValueBox({
valueBox(sum(jabar_merged$JK, na.rm = TRUE),
"Total Kasus TBC",
icon = icon("lungs"),
color = "red")
})
output$max_cases <- renderValueBox({
max_row <- jabar_merged %>%
st_drop_geometry() %>%
filter(JK == max(JK, na.rm = TRUE)) %>%
slice(1)
valueBox(
value = format(max_row$JK, big.mark = ","),
subtitle = paste0("Kasus TBC Tertinggi (", max_row$WADMKK, ")"),
icon = icon("arrow-up"),
color = "red"
)
})
output$num_districts <- renderValueBox({
valueBox(length(unique(data$kab_kota)),
"Jumlah Kabupaten/Kota",
icon = icon("map-marker-alt"),
color = "blue")
})
# Data Table & Download
output$data_table <- renderDT({
data_display <- data
colnames(data_display)[colnames(data_display) == "kab_kota"] <- "Kabupaten/Kota"
colnames(data_display)[colnames(data_display) == "JK"] <- "Jumlah Kasus TBC"
colnames(data_display)[colnames(data_display) == "JP"] <- "Jumlah Penduduk"
datatable(
head(data_display, input$n_rows),
options = list(scrollX = TRUE, pageLength = 10),
rownames = FALSE
)
})
output$download_data <- downloadHandler(
filename = function() {
paste0("Data_TBC_Jabar_", Sys.Date(), ".xlsx")
},
content = function(file) {
write.xlsx(data, file)
}
)
# ============================================================
# 4. STATISTIKA DESKRIPTIF
# ============================================================
# Update pilihan variabel secara dinamis
desc_data <- reactive({
req(input$desc_vars)
data %>%
dplyr::select(all_of(input$desc_vars))
})
desc_data_renamed <- reactive({
desc_data() %>%
rename_with(~ recode(.x,
"JK" = "Jumlah Kasus TBC",
"JP" = "Jumlah Penduduk",
"HIV" = "Kasus HIV",
"DB" = "Penderita Diabetes",
"ST" = "Stunting (%)",
"KP" = "Kepadatan Penduduk",
"Miskin" = "Penduduk Miskin (ribu)",
"Faskes" = "Fasilitas Kesehatan",
"Sanitasi" = "Sanitasi Layak (%)",
"Rate_1k" = "Prevalensi per 1.000",
"RR" = "Relative Risk"
))
})
output$desc_histogram <- renderPlotly({
req(input$desc_vars)
df <- desc_data_renamed()
plots <- lapply(names(df), function(var) {
plot_ly(
df,
x = ~get(var),
type = "histogram",
name = var
) %>%
layout(xaxis = list(title = var),
yaxis = list(title = "Frekuensi"))
})
subplot(plots, nrows = ceiling(length(plots) / 2),
margin = 0.05, shareX = FALSE, shareY = FALSE)
})
output$desc_boxplot <- renderPlotly({
req(input$desc_vars)
df <- desc_data_renamed()
plots <- lapply(names(df), function(var) {
plot_ly(
df,
y = ~get(var),
type = "box",
name = var,
boxmean = TRUE,
marker = list(color = "blue")
) %>%
layout(yaxis = list(title = var))
})
subplot(plots,
nrows = ceiling(length(plots) / 2),
margin = 0.05,
shareX = FALSE,
shareY = FALSE)
})
output$desc_summary <- renderDT({
req(input$desc_vars)
df <- desc_data_renamed()
summary_df <- df %>%
summarise(across(
everything(),
list(
Mean = ~mean(.x, na.rm = TRUE),
SD = ~sd(.x, na.rm = TRUE),
Min = ~min(.x, na.rm = TRUE),
Max = ~max(.x, na.rm = TRUE)
)
)) %>%
pivot_longer(
everything(),
names_to = c("Variabel", "Statistik"),
names_sep = "_",
values_to = "Nilai"
)
datatable(
summary_df,
options = list(pageLength = 10),
rownames = FALSE
)
})
output$download_desc <- downloadHandler(
filename = function() {
paste0("Statistik_Deskriptif_TBC_", Sys.Date(), ".xlsx")
},
content = function(file) {
write.xlsx(desc_data_renamed(), file)
}
)
# ============================================================
# 5. SEBARAN TBC
# ============================================================
output$total_kasus <- renderValueBox({
valueBox(
format(total_JK, big.mark = ","),
"Total Kasus TBC",
icon = icon("procedures"),
color = "red"
)
})
output$max_rr <- renderValueBox({
valueBox(
round(max(jabar_merged$RR, na.rm = TRUE), 2),
"RR Tertinggi",
icon = icon("exclamation-triangle"),
color = "orange"
)
})
output$total_wilayah <- renderValueBox({
valueBox(
nrow(data),
"Kabupaten/Kota",
icon = icon("map"),
color = "blue"
)
})
# PETA SEBARAN JK
output$map_jk <- renderLeaflet({
pal <- colorNumeric("YlOrRd", domain = jabar_merged$JK)
leaflet(jabar_merged) %>%
addTiles() %>%
addPolygons(
fillColor = ~pal(JK),
weight = 1,
opacity = 1,
color = "black",
fillOpacity = 0.7,
label = ~paste0(
WADMKK, "<br>",
"JK: ", round(JK, 2), "<br>",
"Kasus: ", format(JK, big.mark = ",")
),
highlightOptions = highlightOptions(weight = 3, fillOpacity = 0.9)
) %>%
addLegend(pal = pal, values = ~JK,
title = "Jumlah Kasus TBC (JK)",
position = "bottomright")
})
moran_result <- reactive({
x <- jabar_merged$RR
names(x) <- rownames(jabar_sp)
moran.test(x, WL, zero.policy = TRUE)
})
geary_result <- reactive({
x <- jabar_merged$RR
names(x) <- rownames(jabar_sp)
geary.test(x, WL, zero.policy = TRUE)
})
output$moranI_box <- renderValueBox({
m <- moran_result()
valueBox(
round(m$estimate[1], 4),
paste0("Moran's I (p: ", round(m$p.value, 4), ")"),
icon = icon("project-diagram"),
color = ifelse(m$p.value < 0.05, "green", "red")
)
})
output$gearyC_box <- renderValueBox({
g <- geary_result()
valueBox(
round(g$estimate[1], 4),
paste0("Geary's C (p: ", round(g$p.value, 4), ")"),
icon = icon("chart-area"),
color = ifelse(g$p.value < 0.05, "green", "red")
)
})
df_moran <- reactive({
req(model_ols())
res <- residuals(model_ols())
data.frame(
x = as.numeric(scale(res)),
wx = as.numeric(scale(lag.listw(WL, res, zero.policy = TRUE)))
)
})
output$moran_plot <- renderPlotly({
df <- df_moran()
fit <- lm(wx ~ x, data = df)
colors <- c("Q1: High-High" = "#d62728",
"Q2: Low-High" = "#ff7f0e",
"Q3: Low-Low" = "#1f77b4",
"Q4: High-Low" = "#2ca02c")
plot_ly(df, x = ~x, y = ~wx) %>%
add_markers() %>%
add_lines(x = ~x, y = fitted(fit)) %>%
plotly::layout(
title = "Moran's I Scatterplot",
xaxis = list(title = "Residual"),
yaxis = list(title = "Spatial Lag Residual")
)
})
# LISA Map - Versi Perbaikan
output$lisa_map <- renderLeaflet({
# 1. Ambil data RR dan pastikan numerik
x <- as.numeric(jabar_merged$RR)
# 2. Hitung Local Moran
# WL didapat dari global section (listw)
lisa <- localmoran(x, WL, zero.policy = TRUE)
# 3. Scaling untuk klasifikasi kuadran
x_std <- as.numeric(scale(x))
x_lag_std <- as.numeric(scale(lag.listw(WL, x, zero.policy = TRUE)))
# 4. Inisialisasi Kategori
lisa_class <- rep("Not Significant", length(x))
p_val <- lisa[, 5] # Kolom p-value
# 5. Klasifikasi berdasarkan Signifikansi (p < 0.05)
sig <- p_val < 0.05
lisa_class[sig & x_std > 0 & x_lag_std > 0] <- "High-High"
lisa_class[sig & x_std < 0 & x_lag_std < 0] <- "Low-Low"
lisa_class[sig & x_std > 0 & x_lag_std < 0] <- "High-Low"
lisa_class[sig & x_std < 0 & x_lag_std > 0] <- "Low-High"
# 6. Masukkan kembali ke data frame (PENTING: Gunakan faktor agar legenda konsisten)
lisa_levels <- c("High-High", "Low-Low", "High-Low", "Low-High", "Not Significant")
jabar_merged$lisa_cat <- factor(lisa_class, levels = lisa_levels)
# 7. Definisi Palet Warna sesuai standar LISA
lisa_colors <- c("#e31a1c", "#1f78b4", "#fb9a99", "#a6cee3", "#eeeeee") # Red, Blue, Pink, LightBlue, Grey
pal <- colorFactor(palette = lisa_colors, domain = jabar_merged$lisa_cat, levels = lisa_levels)
# 8. Render Peta
leaflet(jabar_merged) %>%
addTiles() %>%
addPolygons(
fillColor = ~pal(lisa_cat),
weight = 1,
opacity = 1,
color = "black",
fillOpacity = 0.8,
label = ~paste0(WADMKK, ": ", lisa_cat, " (p-val: ", round(p_val, 3), ")"),
highlightOptions = highlightOptions(weight = 3, color = "white", bringToFront = TRUE)
) %>%
addLegend(
pal = pal,
values = ~lisa_cat,
title = "LISA Cluster",
position = "bottomright"
)
})
# Getis-Ord Gi*
output$getisord_map <- renderLeaflet({
gi <- localG(jabar_merged$RR, WL, zero.policy = TRUE)
jabar_merged$z_Gistar <- as.numeric(gi)
jabar_merged$hotcold <- ifelse(
jabar_merged$z_Gistar > 1.96, "Hot Spot",
ifelse(jabar_merged$z_Gistar < -1.96, "Cold Spot", "Not Significant")
)
pal <- colorFactor(
palette = c(
"Hot Spot" = "#b2182b",
"Cold Spot" = "#2166ac",
"Not Significant" = "grey85"
),
domain = jabar_merged$hotcold
)
leaflet(jabar_merged) %>%
addTiles() %>%
addPolygons(
fillColor = ~pal(hotcold),
weight = 1,
opacity = 1,
color = "black",
fillOpacity = 0.7,
label = ~paste0(
WADMKK, ": ", hotcold,
"<br>z-score Gi*: ", round(z_Gistar, 2),
"<br>RR: ", round(RR, 2),
"<br>Kasus TBC: ", format(JK, big.mark = ",")
),
highlightOptions = highlightOptions(weight = 3, fillOpacity = 0.9)
) %>%
addLegend(
pal = pal,
values = ~hotcold,
title = "Getis-Ord Gi* (RR)",
position = "bottomright"
)
})
# ============================================================
# 6. UKURAN EPIDEMIOLOGI
# ============================================================
output$prev_jabar_box <- renderValueBox({
valueBox(
value = round(prev_jabar_1k, 2),
subtitle = "Prevalensi TBC Jawa Barat (per 1.000)",
icon = icon("chart-line"),
color = "light-blue"
)
})
output$rr_jabar_box <- renderValueBox({
valueBox(
value = RR_jabar,
subtitle = "Relative Risk (Baseline Provinsi)",
icon = icon("balance-scale"),
color = "green"
)
})
# PETA PREVALENSI TBC (per 1.000 penduduk)
output$map_prevalensi <- renderLeaflet({
breaks_prev <- c(0, 2, 4, 6, Inf)
labels_prev <- c("<2", "2-4", "4-6", ">6")
jabar_merged$Kategori_Prev <- cut(
jabar_merged$Prev_1k,
breaks = breaks_prev,
labels = labels_prev,
right = TRUE,
include.lowest = TRUE
)
pal <- colorFactor(
palette = c("< 2" = "#ffffb2",
"2-4" = "#fecc5c",
"4-6" = "#fd8d3c",
"> 6" = "#e31a1c"),
domain = jabar_merged$Kategori_Prev
)
leaflet(jabar_merged) %>%
addTiles() %>%
addPolygons(
fillColor = ~pal(Kategori_Prev),
weight = 1,
opacity = 1,
color = "black",
fillOpacity = 0.7,
label = ~paste0(
WADMKK,
": ", round(Prev_1k, 2),
" per 1.000 penduduk"
),
highlightOptions = highlightOptions(weight = 3, fillOpacity = 0.9)
) %>%
addLegend(
pal = pal,
values = ~Kategori_Prev,
title = "Prevalensi TBC (per 1.000)",
position = "bottomright"
)
})
# Prevalensi per Wilayah (Tabel & Barplot)
prev_wilayah <- reactive({
req(input$prev_cities)
jabar_merged %>%
filter(WADMKK %in% input$prev_cities) %>%
# Gunakan dplyr:: sebelum select untuk menghindari konflik
dplyr::select(WADMKK, JK, JP, Prev_1k) %>%
st_drop_geometry() %>%
arrange(desc(Prev_1k))
})
output$prev_table <- renderDT({
df <- prev_wilayah()
datatable(
df,
colnames = c(
"Kabupaten/Kota",
"Jumlah Kasus TBC",
"Jumlah Penduduk",
"Prevalensi (per 1.000)"
),
options = list(dom = 't', paging = FALSE, ordering = FALSE)
) %>%
formatRound(columns = "Prev_1k", digits = 2) %>%
formatCurrency(columns = c("JK", "JP"), currency = "", digits = 0)
})
output$prev_bar <- renderPlotly({
df <- prev_wilayah()
plot_ly(
df,
x = ~Prev_1k,
y = ~reorder(WADMKK, Prev_1k),
type = 'bar',
orientation = 'h',
marker = list(
color = '#3498db',
line = list(color = 'white', width = 1)
),
text = ~round(Prev_1k, 2),
textposition = 'outside',
hoverinfo = 'text',
hovertext = ~paste0(
"<b>", WADMKK, "</b><br>",
"Prevalensi: ", round(Prev_1k, 2), " per 1.000<br>",
"Kasus: ", format(JK, big.mark = ","), "<br>",
"Penduduk: ", format(JP, big.mark = ",")
)
) %>%
plotly::layout(
title = "Prevalensi TBC per Kabupaten/Kota",
xaxis = list(title = "Prevalensi (per 1.000 penduduk)"),
yaxis = list(title = ""),
margin = list(l = 150)
)
})
# Relative Risk (RR) per Wilayah
rr_wilayah <- reactive({
req(input$prev_cities)
jabar_merged %>%
filter(WADMKK %in% input$prev_cities) %>%
# Gunakan dplyr:: sebelum select
dplyr::select(WADMKK, JK, JP, RR) %>%
st_drop_geometry() %>%
arrange(desc(RR))
})
output$rr_table <- renderDT({
df <- rr_wilayah()
datatable(
df,
colnames = c(
"Kabupaten/Kota",
"Jumlah Kasus TBC",
"Jumlah Penduduk",
"Relative Risk (RR)"
),
options = list(dom = 't', paging = FALSE, ordering = FALSE)
) %>%
formatRound(columns = "RR", digits = 2) %>%
formatCurrency(columns = c("JK", "JP"), currency = "", digits = 0)
})
# PETA RR
output$map_rr <- renderLeaflet({
# Klasifikasi RR (interpretatif & umum dipakai)
breaks_rr <- c(0, 0.75, 1.00, 1.25, Inf)
labels_rr <- c("< 0.75 (Risiko Rendah)",
"0.75-1.00 (= Rata-rata Provinsi)",
"1.01-1.25 (Risiko Tinggi)",
"> 1.25 (Risiko Sangat Tinggi)")
jabar_merged$Kategori_RR <- cut(
jabar_merged$RR,
breaks = breaks_rr,
labels = labels_rr,
right = TRUE,
include.lowest = TRUE
)
pal <- colorFactor(
palette = c("< 0.75 (Risiko Rendah)" = "#ffffb2",
"0.75 – 1.00 (≈ Rata-rata)" = "#fecc5c",
"1.01 – 1.25 (Risiko Tinggi)" = "#fd8d3c",
"> 1.25 (Risiko Sangat Tinggi)" = "#e31a1c"),
domain = jabar_merged$Kategori_RR
)
leaflet(jabar_merged) %>%
addTiles() %>%
addPolygons(
fillColor = ~pal(Kategori_RR),
weight = 1,
opacity = 1,
color = "black",
fillOpacity = 0.7,
label = ~paste0(
WADMKK, ": RR = ", round(RR, 2)
),
highlightOptions = highlightOptions(
weight = 3,
fillOpacity = 0.9
)
) %>%
addLegend(
pal = pal,
values = ~Kategori_RR,
title = "Relative Risk (RR)",
position = "bottomright"
)
})
# RR per Wilayah
rr_wilayah <- reactive({
req(input$prev_cities)
jabar_merged %>%
filter(WADMKK %in% input$prev_cities) %>%
select(WADMKK, JK, JP, RR) %>%
st_drop_geometry() %>%
arrange(desc(RR))
})
output$rr_table <- renderDT({
df <- rr_wilayah()
datatable(
df,
colnames = c(
"Kabupaten/Kota",
"Jumlah Kasus TBC",
"Jumlah Penduduk",
"Relative Risk (RR)"
),
options = list(dom = 't', paging = FALSE, ordering = FALSE)
) %>%
formatRound(columns = "RR", digits = 2) %>%
formatCurrency(columns = c("JK", "JP"), currency = "", digits = 0)
})
# BARPLOT RR PER WILAYAH
output$rr_bar <- renderPlotly({
req(input$prev_cities)
df <- rr_wilayah()
plot_ly(
df,
x = ~RR,
y = ~reorder(WADMKK, RR),
type = 'bar',
orientation = 'h',
mode = "markers",
marker = list(
color = '#2ecc71', # hijau → risiko
line = list(color = 'white', width = 1)
),
text = ~round(RR, 2),
textposition = 'outside',
hoverinfo = 'text',
hovertext = ~paste0(
"<b>", WADMKK, "</b><br>",
"RR: ", round(RR, 2), "<br>",
"Kasus TBC: ", format(JK, big.mark = ","), "<br>",
"Penduduk: ", format(JP, big.mark = ",")
)
) %>%
plotly::layout(
title = "Relative Risk (RR) TBC per Kabupaten/Kota",
xaxis = list(title = "Relative Risk (RR)"),
yaxis = list(title = ""),
margin = list(l = 150)
)
})
# ============================================================
# Pemodelan
# ============================================================
# --- 7.1 REGRESI OLS ---
model_ols <- reactive({
lm(Rate_1k ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
data = jabar_merged)
})
output$ols_summary <- renderPrint({
summary(model_ols())
})
output$ols_r2 <- renderValueBox({
model_lm <- lm(Rate_1k ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP, data = jabar_merged)
valueBox(round(summary(model_lm)$r.squared, 4), "R-Squared", icon = icon("chart-line"), color = "blue")
})
output$ols_rmse <- renderValueBox({
model_lm <- lm(Rate_1k ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP, data = jabar_merged)
rmse_val <- sqrt(mean(residuals(model_lm)^2))
valueBox(round(rmse_val, 4), "RMSE", icon = icon("calculator"), color = "red")
})
output$ols_asumsi_table <- renderDT({
model_lm <- lm(Rate_1k ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP, data = jabar_merged)
# Dataframe hasil uji asumsi dari syntax Anda
asumsi_df <- data.frame(
Uji = c("Normalitas (Jarque-Bera)", "Homoskedastisitas (BP-Test)", "Multikolinearitas (Max VIF)"),
Statistik = c(round(tseries::jarque.bera.test(residuals(model_lm))$statistic, 3),
round(lmtest::bptest(model_lm)$statistic, 3),
round(max(car::vif(model_lm)), 3)),
P_Value = c(round(tseries::jarque.bera.test(residuals(model_lm))$p.value, 4),
round(lmtest::bptest(model_lm)$p.value, 4),
"-")
)
datatable(asumsi_df, options = list(dom = 't'), rownames = FALSE)
})
# --- 7.2 POISSON & NEGATIVE BINOMIAL ---
output$irr_pg_table <- renderDT({
table6 <- data.frame(
Variabel = c("HIV", "Stunting", "Diabetes", "Sanitasi", "Faskes", "Miskin", "Kepadatan"),
Nilai_IRR = c(1.0006, 1.0420, 0.99997, 1.0146, 1.0046, 1.0030, 1.00005),
Interpretasi = c(
"Tidak signifikan. Tidak memengaruhi laju TBC setelah koreksi overdispersion.",
"Tidak signifikan. Arah efek positif, namun tidak cukup bukti statistik.",
"Signifikan. Setiap tambahan diabetes menurunkan laju TBC sangat kecil.",
"Tidak signifikan. Tidak berpengaruh setelah heterogenitas diperhitungkan.",
"Tidak signifikan. Efek tidak kuat setelah koreksi dispersi.",
"Tidak signifikan. Efek sosial melemah setelah overdispersion ditangani.",
"Tidak signifikan. Pengaruh kepadatan tidak konsisten."
)
)
datatable(table6, rownames = FALSE, options = list(dom = 't', scrollX = TRUE)) %>%
formatRound(columns = "Nilai_IRR", digits = 5)
})
# --- 7.3 BYM-INLA ---
output$bym_ci_table <- renderDT({
table7 <- data.frame(
Variabel = c("Intercept", "HIV", "Stunting", "Diabetes", "Sanitasi", "Faskes", "Miskin", "Kepadatan"),
Mean = c(-1.066, 0.000, 0.021, 0.000, 0.007, 0.005, 0.002, 0.000),
CI_IRR = c("-3.500 - 1.367", "-0.002 - 0.003", "-0.058 - 0.100", "0.000 - 0.000",
"-0.016 - 0.030", "-0.027 - 0.037", "-0.006 - 0.009", "0.000 - 0.000"),
Interpretasi = c(
"Nilai dasar log-risiko TBC. Tidak signifikan karena interval mencakup nol.",
"Tidak ada bukti kuat HIV berpengaruh setelah mempertimbangkan efek spasial.",
"Stunting tidak menunjukkan pengaruh signifikan terhadap risiko TBC.",
"Tidak terdapat pengaruh bermakna diabetes dalam model spasial.",
"Akses sanitasi tidak berhubungan signifikan setelah efek spasial dihitung.",
"Faskes tidak berpengaruh signifikan terhadap TBC antar wilayah.",
"Jumlah penduduk miskin tidak menunjukkan hubungan signifikan.",
"Kepadatan penduduk tidak berpengaruh signifikan dalam model BYM."
)
)
datatable(table7, rownames = FALSE, options = list(dom = 't', scrollX = TRUE)) %>%
formatRound(columns = "Mean", digits = 3)
})
output$map_bym_rr <- renderLeaflet({
# model_bym harus sudah di-run di server/global
# RR dihitung dari: model_bym$summary.fitted.values$mean / jabar_merged$E
pal <- colorNumeric("YlOrRd", domain = jabar_merged$RR)
leaflet(jabar_merged) %>%
addTiles() %>%
addPolygons(
fillColor = ~pal(RR),
weight = 1, color = "black", fillOpacity = 0.7,
label = ~paste0(WADMKK, ": RR Bayesian = ", round(RR, 3))
) %>%
addLegend(pal = pal, values = ~RR, title = "RR (BYM-INLA)")
})
output$bym_dic <- renderValueBox({
# Asumsi model_bym sudah tersedia
valueBox(round(350.42, 2), "DIC (BYM Model)", icon = icon("info-circle"), color = "purple")
})
# ============================================================
# Analisis Spasial
# ============================================================
# OLS
output$ols_formula <- renderUI({
withMathJax("
$$ JK_i = 2368 + 0.0005013 JP_i + 2.195 HIV_i - 181.9 ST_i
- 0.08767 DB_i - 48.94 Sanitasi_i + 139.9 Faskes_i
+ 21.98 Miskin_i + 0.2768 KP_i + \\varepsilon_i $$
")
})
output$ols_table <- renderTable({
req(model_ols())
broom::tidy(model_ols())
})
output$moran_ols <- renderTable({
df_moran_res <- df_moran()
cor(df_moran_res$x, df_moran_res$wx) # contoh sederhana
})
output$ols_interpretation <- renderText({
"Model OLS memiliki R-square sebesar 0,912 yang menunjukkan bahwa 91,2% variasi jumlah kasus TBC
dapat dijelaskan oleh variabel independen dalam model. Uji Moran's I pada residual
menunjukkan tidak adanya autokorelasi spasial yang signifikan."
})
#SAR
model_spatial <- reactive({
lagsarlm(
Rate_1k ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
data = jabar_sp,
listw = WL,
zero.policy = TRUE
)
})
output$spatial_classic_table <- renderPrint({
summary(model_spatial())
})
model_ols_sp <- reactive({
lm(
Rate_1k ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
data = jabar_sp@data
)
})
lm_test_result <- reactive({
lm.RStests(
model_ols_sp(),
WL,
test = c("LMlag", "LMerr", "RLMlag", "RLMerr"),
zero.policy = TRUE
)
})
output$lm_test <- renderDT({
rst <- lm_test_result()
datatable(
data.frame(
Test = names(rst),
Statistic = sapply(rst, function(x) x$statistic),
DF = sapply(rst, function(x) x$parameter),
P_value = sapply(rst, function(x) x$p.value)
),
options = list(dom = "t"),
rownames = FALSE
)
})
output$spatial_classic_interp <- renderText({
"Nilai p-value untuk Rho lebih besar dari 0.05 (tidak signifikan). Artinya tidak terdapat ketergantungan spasial (autokorelasi spasial) yang signifikan pada variabel respon (Rate_1k). Dengan kata lain, tingkat kejadian di suatu wilayah tidak dipengaruhi secara signifikan oleh tingkat kejadian di wilayah tetangganya.
Hasil Uji diagnostik (LM Test) menunjukkan tidak ada dependensi spasial"
})
# INTERPOLASI
output$interp_table <- renderTable({
interp_summary # ringkasan IDW, Trend Surface, Kriging
})
output$interp_interp <- renderText({
"Metode interpolasi spasial digunakan untuk mengestimasi nilai TBC pada lokasi
yang tidak teramati. Ordinary Kriging memberikan hasil paling stabil karena
mempertimbangkan struktur variogram spasial."
})
#PLOT IDW
output$idw_plot <- renderPlot({
plot(
idw_raster,
main = "Interpolasi IDW Prevalensi TBC",
col = viridis::viridis(100)
)
plot(st_geometry(jabar_merged), add = TRUE)
points(centroid_sp, pch = 20, cex = 0.7)
})
output$trend_plot <- renderPlot({
plot(
trend_raster,
main = "Trend Surface Prevalensi TBC",
col = viridis::viridis(100)
)
plot(st_geometry(jabar_merged), add = TRUE)
})
#KRIGING
output$kriging_plot <- renderPlot({
plot(
kriging_raster,
main = "Ordinary Kriging Prevalensi TBC",
col = viridis::viridis(100)
)
plot(st_geometry(jabar_merged), add = TRUE)
})
output$idw_plot <- renderPlot({
plot(idw_raster, main = "Interpolasi IDW (Prevalensi TBC)")
plot(jabar_merged$geometry, add = TRUE)
})
output$trend_plot <- renderPlot({
plot(trend_raster, main = "Trend Surface (Prevalensi TBC)")
plot(jabar_merged$geometry, add = TRUE)
})
output$kriging_plot <- renderPlot({
plot(kriging_raster, main = "Ordinary Kriging (Prevalensi TBC)")
plot(jabar_merged$geometry, add = TRUE)
})
}
# ============================================================
# 4. RUN APP
# ============================================================
shinyApp(ui = ui, server = server)