Syntax Dashboard Epidemiologi Tuberkulosis 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

---

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(performance)
library(dplyr)

# ============================================================
# 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
data <- data %>%
  mutate(
    JK = as.numeric(JK),
    JP = as.numeric(JP)
  )

# Hitung risk rata-rata Jawa Barat
risk_jabar <- sum(data$JK, na.rm = TRUE) /
  sum(data$JP, na.rm = TRUE)

# Gabungkan & hitung RR
jabar_merged <- Indo %>%
  left_join(data, by = c("WADMKK" = "kab_kota")) %>%
  mutate(
    Prev_1k = (JK / JP) * 1000,
    RR = (JK / JP) / risk_jabar
  )

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

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"))
    )
  ),
  
  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")
                  )
                )
              )
            )
    ),
    # =====================================================
    # 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")
                )
              )
            )
    )
    
  )
)


# ============================================================
# 3. SERVER SECTION
# ============================================================

server <- function(input, output, session) {
  
  # ============================================================
  # 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({
    df <- desc_data_renamed()
    
    numeric_vars <- df %>%
      dplyr::select(where(is.numeric)) %>%
      colnames()
    
    updatePickerInput(
      session,
      "desc_vars",
      choices = numeric_vars,
      selected = numeric_vars
    )
  })
  
  # ============================================================
  # 2. DATA OVERVIEW
  # ============================================================
  # ValueBoxes
  output$total_cases <- renderValueBox({
    valueBox(sum(data$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
  # ============================================================
  
  desc_data <- reactive({
    data   # langsung pakai seluruh data tanpa filtering
  })
  
  desc_data_renamed <- reactive({
    desc_data() %>%
      rename(
        "Kabupaten/Kota" = kab_kota,
        "Jumlah Kasus TBC" = JK,
        "Jumlah Penduduk" = JP,
        "Kasus HIV" = HIV,
        "Penderita Diabetes" = DB,
        "Stunting (%)" = ST,
        "Kepadatan Penduduk" = KP,
        "Penduduk Miskin (ribu)" = Miskin,
        "Fasilitas Kesehatan" = Faskes,
        "Sanitasi Layak (%)" = Sanitasi,
        "Prevalensi per 1.000" = Rate_1k,
        "Relative Risk" = RR
      )
  })
  output$desc_histogram <- renderPlotly({
    req(input$desc_vars)
    df <- desc_data_renamed()
    
    plots <- lapply(input$desc_vars, function(var) {
      plot_ly(
        df,
        x = ~get(var),
        type = "histogram",
        name = var) %>%
        plotly::layout(
          xaxis = list(title = var),
          yaxis = list(title = "Frekuensi")
        )
    })
    
    subplot(plots, nrows = ceiling(length(input$desc_vars) / 2), 
            margin = 0.05, shareX = FALSE, shareY = FALSE)
  })
  
  output$desc_boxplot <- renderPlotly({
    req(input$desc_vars)
    df <- desc_data_renamed()
    
    plots <- lapply(input$desc_vars, function(var) {
      plot_ly(
        df,
        y = ~get(var),
        type = "box",
        name = var,
        boxmean = TRUE,
        marker = list(color = "blue")
      ) %>%
        plotly::layout(yaxis = list(title = var))
    })
    
    subplot(
      plots,
      nrows = ceiling(length(input$desc_vars) / 2),
      margin = 0.05,
      shareX = FALSE,
      shareY = FALSE
    )
  })
  output$desc_summary <- renderDT({
    req(input$desc_vars)
    
    # 1. Ambil data dasar
    df <- desc_data_renamed()
    
    # 2. Lakukan perhitungan statistik
    # Gunakan dplyr::select untuk menghindari bentrok package lain
    summary_df <- df %>%
      select(all_of(input$desc_vars)) %>%
      summarise(across(everything(), list(
        Min = ~min(., na.rm = TRUE),
        Q1 = ~quantile(., 0.25, na.rm = TRUE),
        Median = ~median(., na.rm = TRUE),
        Mean = ~mean(., na.rm = TRUE),
        Q3 = ~quantile(., 0.75, na.rm = TRUE),
        Max = ~max(., na.rm = TRUE),
        SD = ~sd(., na.rm = TRUE)
      ), .names = "{.col}###{.fn}")) %>% # Gunakan separator unik
      pivot_longer(
        everything(),
        names_to = "Variable_Stat",
        values_to = "Value"
      ) %>%
      separate(Variable_Stat, into = c("Variable", "Statistic"), sep = "###") %>%
      pivot_wider(names_from = Statistic, values_from = Value)
    
    # 3. Render ke Tabel
    datatable(
      summary_df,
      rownames = FALSE,
      options = list(
        dom = "t", 
        paging = FALSE, 
        scrollX = TRUE,
        columnDefs = list(list(className = 'dt-center', targets = "_all"))
      )
    ) %>%
      formatRound(columns = 2:8, digits = 2)
  })
  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")
    )
  })
  
  output$moran_plot <- renderPlotly({
    req(jabar_merged$RR)
    x <- as.numeric(scale(jabar_merged$RR))
    wx <- as.numeric(scale(lag.listw(WL, jabar_merged$RR, zero.policy = TRUE)))
    
    df_moran <- data.frame(x = x, wx = wx, region = jabar_merged$WADMKK)
    fit <- lm(wx ~ x, data = df_moran)
    
    
    colors <- c("Q1: High-High" = "#d62728",
                "Q2: Low-High" = "#ff7f0e",
                "Q3: Low-Low" = "#1f77b4",
                "Q4: High-Low" = "#2ca02c")
    
    plot_ly(df_moran, x = ~x, y = ~wx) %>%
      add_markers(text = ~region, hoverinfo = "text", name = "Wilayah",
                  marker = list(color = '#1f77b4', size = 8)) %>%
      add_lines(x = ~x, y = fitted(fit), name = "Regresi",
                line = list(color = "red", width = 2)) %>%
      plotly::layout(title = "Moran's I Scatterplot",
                     xaxis = list(title = "Standardized RR (x)"),
                     yaxis = list(title = "Spatial Lag RR (Wx)"),
                     showlegend = FALSE)
    
  })
  
  # 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',
      mode = "markers",
      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 ---
  output$ols_summary <- renderPrint({
    model_lm <- lm(Rate_1k ~ HIV + ST + DB + Sanitasi + Faskes + Miskin + KP, data = jabar_merged)
    summary(model_lm)
  })
  
  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") 
  })
}

# ============================================================
# 4. RUN APP
# ============================================================

shinyApp(ui = ui, server = server)