Backlog Kelayakhunian 2025

knitr::opts_chunk$set(echo = TRUE)

Data yang digunakan adalah Susenas Maret 2025.

1 Backlog Kelayakhunian

Backlog kelayakhunian merupakan indikator yang menggambarkan kondisi rumah tangga yang menempati bangunan tempat tinggal milik sendiri, namun dalam kondisi yang tidak layak huni.

2 Set Design Sampling

Design sampling yang digunakan mengikuti design sampling yang dilakukan pada SPSS yaitu Cluster Sampling design with replacement.

susenas.design<- svydesign(id=~psu, strata=~strata, data = dataku, weights=~fwt)
susenas.design

## Stratified 1 - level Cluster Sampling design (with replacement)
## With (34473) clusters.
## svydesign(id = ~psu, strata = ~strata, data = dataku, weights = ~fwt)

summary(susenas.design$prob)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000201 0.003795 0.008842 0.029646 0.025600 0.968267

3 Set Output Statistics

Kualitas presisi hasil estimasi suatu survei bisa diamati dari nilai RSE yang dihasilkan. Kesalahan sampling dari beberapa estimasi harus digunakan secara hati-hati. Untuk estimasi yang berdasarkan jumlah kasus yang kecil, kesalahan relatif cenderung besar.

Batasan nilai RSE ( Australian Bureau Statistics):

1. Jika \(RSE \le 25\%\), estimasi bersifat presisi.

2. Jika \(25\% < RSE \le 50\%\), estimasi perlu dilakukan dengan hati-hati.

3. Jika \(RSE > 50\%\), estimasi dianggap sangat tidak presisi.

Provinsi

options(survey.adjust.domain.lonely=TRUE)
options(survey.lonely.psu="adjust")
hasilP = svyby(formula = ~BACK_LOG, denom= ~denom, ~prov, design = susenas.design, data = dataku, deff=TRUE, svyratio, vartype=c("se","ci","ci","cv","cvpct","var"))
hasilP[is.na(hasilP)] <- 0
hasilP$theta = round(hasilP$`BACK_LOG/denom`*100,2)
hasilP$SE = round(hasilP$`se.BACK_LOG/denom`*100,2)
hasilP$VAR = round(hasilP$SE*hasilP$SE,2)
hasilP$CI_LOWER = round(hasilP$`ci_l`*100,2)
hasilP$CI_LOWER[hasilP$CI_LOWER<0] <- 0
hasilP$CI_UPPER = round(hasilP$`ci_u`*100,2)
hasilP$RSE = round(hasilP$`cv%`,2)
hasilP$DEFF = round(hasilP$DEff,2)

Kabupaten/kota

options(survey.adjust.domain.lonely=TRUE)
options(survey.lonely.psu="adjust")
hasil = svyby(formula = ~BACK_LOG, denom= ~denom, ~kabu, design = susenas.design, data = dataku, deff=TRUE, svyratio, vartype=c("se","ci","ci","cv","cvpct","var"))
hasil[is.na(hasil)] <- 0
hasil$theta = round(hasil$`BACK_LOG/denom`*100,2)
hasil$SE = round(hasil$`se.BACK_LOG/denom`*100,2)
hasil$VAR = round(hasil$SE*hasil$SE,2)
hasil$CI_LOWER = round(hasil$`ci_l`*100,2)
hasil$CI_LOWER[hasil$CI_LOWER<0] <- 0
hasil$CI_UPPER = round(hasil$`ci_u`*100,2)
hasil$RSE = round(hasil$`cv%`,2)
hasil$DEFF = round(hasil$DEff,2)

4 Level Estimate

4.1 Estimation by province

outputP = as.data.frame(cbind(hasilP$prov, hasilP$theta, hasilP$SE, hasilP$VAR, hasilP$CI_LOWER, hasilP$CI_UPPER, hasilP$RSE, hasilP$DEFF))
names(outputP) = c("Prov","Estimasi","SE","VAR","CI LOWER","CI UPPER","RSE","DEFF")
(outputP[order(outputP$Prov, decreasing = FALSE), ]   )

##    Prov Estimasi   SE  VAR CI LOWER CI UPPER  RSE DEFF
## 1    11    24.95 0.66 0.44    23.66    26.25 2.64 3.21
## 2    12    16.02 0.56 0.31    14.93    17.11 3.46 4.92
## 3    13    26.61 0.85 0.72    24.95    28.28 3.19 4.32
## 4    14    15.37 0.76 0.58    13.89    16.86 4.93 3.80
## 5    15    28.16 1.18 1.39    25.85    30.46 4.18 4.96
## 6    16    28.16 0.98 0.96    26.24    30.08 3.48 5.64
## 7    17    37.37 1.27 1.61    34.89    39.85 3.39 4.12
## 8    18    28.00 1.04 1.08    25.97    30.03 3.71 5.78
## 9    19    59.86 1.37 1.88    57.17    62.55 2.29 3.36
## 10   21    31.00 2.16 4.67    26.77    35.23 6.97 9.11
## 11   31    28.05 1.12 1.25    25.85    30.25 4.00 3.18
## 12   32    33.00 0.69 0.48    31.66    34.35 2.08 5.27
## 13   33    22.63 0.51 0.26    21.62    23.63 2.27 4.41
## 14   34     8.24 0.82 0.67     6.64     9.84 9.91 3.48
## 15   35    20.17 0.49 0.24    19.22    21.12 2.41 4.61
## 16   36    23.94 1.08 1.17    21.82    26.07 4.53 4.57
## 17   51     8.93 0.79 0.62     7.39    10.48 8.82 4.79
## 18   52    25.43 1.03 1.06    23.40    27.45 4.06 4.01
## 19   53    34.56 0.87 0.76    32.86    36.26 2.51 4.38
## 20   61    24.37 1.03 1.06    22.35    26.40 4.24 5.43
## 21   62    28.94 1.18 1.39    26.63    31.26 4.08 5.67
## 22   63    30.41 1.07 1.14    28.32    32.50 3.50 4.69
## 23   64    16.90 1.19 1.42    14.57    19.23 7.03 6.55
## 24   65    19.09 1.57 2.46    16.00    22.17 8.25 4.31
## 25   71    18.99 0.87 0.76    17.29    20.70 4.57 4.22
## 26   72    30.81 1.02 1.04    28.81    32.82 3.32 3.88
## 27   73    17.88 0.63 0.40    16.65    19.12 3.52 4.31
## 28   74    17.39 0.74 0.55    15.94    18.84 4.26 3.61
## 29   75    20.28 1.15 1.32    18.02    22.54 5.69 2.97
## 30   76    31.11 1.55 2.40    28.06    34.15 4.99 4.07
## 31   81    25.83 1.09 1.19    23.70    27.97 4.21 3.76
## 32   82    23.41 1.19 1.42    21.07    25.75 5.10 4.20
## 33   91    34.72 2.04 4.16    30.72    38.72 5.88 6.28
## 34   92    36.35 2.18 4.75    32.06    40.63 6.01 6.07
## 35   94    32.51 1.75 3.06    29.09    35.94 5.38 6.18
## 36   95    44.85 2.30 5.29    40.34    49.36 5.13 4.43
## 37   96    66.93 1.43 2.04    64.12    69.74 2.14 4.07
## 38   97    89.14 1.32 1.74    86.55    91.72 1.48 7.47

Kelompok RSE Provinsi

hasilP <- hasilP %>%
  mutate(
    RSE = if_else(RSE == 0, NA_real_, RSE)
  ) %>%
  arrange(RSE) %>%
  mutate(
    Urut = row_number(),
    RSEP_kat = case_when(
      is.na(RSE)                 ~ NA_character_,
      RSE < 25.99                ~ "< 25%",
      RSE <= 50.99               ~ "26–50%",
      RSE > 50                   ~ "> 50%"
    )
  )

table(hasilP$RSEP_kat, useNA = "ifany")

## 
## < 25% 
##    38

Visualisasi Provinsi

ggplot(hasilP, aes(x = Urut, y = RSE, group = 1, color = RSEP_kat)) +
  geom_point(size = 3) +
  
  geom_hline(yintercept = 25, linetype = "dashed", color = "red") +
  geom_hline(yintercept = 50, linetype = "dashed", color = "blue") +
  
scale_color_manual(
  name = "Kategori RSE",
  limits = c("< 25%", "26–50%", "> 50%"),
  values = c(
    "< 25%"  = "#1B9E77",
    "26–50%" = "#FF8C00",
    "> 50%"  = "#7570B3"
  ),
  na.value = "#D73027",
  drop = FALSE
) +
    labs(
    title = "RSE BACK_LOG per Provinsi",
    subtitle = "dengan garis batas RSE 25% dan 50%",
    x = "Provinsi",
    y = "Nilai RSE (%)"
  ) +
  
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    legend.position = "bottom",
    axis.text.x = element_text(angle = 40, vjust = 1, hjust = 1)
  )

4.2 Estimation by district

output = as.data.frame(cbind(hasil$kabu, hasil$theta, hasil$SE, hasil$VAR, hasil$CI_LOWER, hasil$CI_UPPER, hasil$RSE, hasil$DEFF))
names(output) = c("Kako","Estimasi","SE","VAR","CI LOWER","CI UPPER","RSE","DEFF")
head((output[order(output$Kako, decreasing = FALSE), ]   ))

##   Kako Estimasi   SE   VAR CI LOWER CI UPPER   RSE DEFF
## 1 1101    34.28 3.20 10.24    28.02    40.55  9.32 2.35
## 2 1102    44.92 4.40 19.36    36.29    53.55  9.80 4.18
## 3 1103    26.78 2.87  8.24    21.15    32.41 10.72 2.63
## 4 1104    30.61 3.63 13.18    23.50    37.72 11.86 3.81
## 5 1105    36.75 3.25 10.56    30.38    43.12  8.85 3.15
## 6 1106    16.34 2.35  5.52    11.73    20.96 14.39 2.43

Kelompok RSE Kabupaten/kota

hasil <- hasil %>%
  mutate(
    RSE = if_else(RSE == 0, NA_real_, RSE)
  ) %>%
  arrange(RSE) %>%
  mutate(
    Urut = row_number(),
    RSE_kat = case_when(
      is.na(RSE)                 ~ NA_character_,
      RSE < 25.99                ~ "< 25%",
      RSE <= 50.99               ~ "26–50%",
      RSE > 50                   ~ "> 50%"
    )
  )

table(hasil$RSE_kat, useNA = "ifany")

## 
##  < 25%  > 50% 26–50% 
##    495      1     18

hasil_plot <- hasil %>%
  mutate(
    RSE_plot = if_else(is.na(RSE), -2, RSE)  # taruh NA di bawah sumbu
  )

ggplot(hasil_plot, aes(x = Urut, y = RSE_plot, color = RSE_kat)) +
  geom_jitter(width = 0.35, height = 0, size = 1.8, alpha = 0.85) +

  geom_hline(yintercept = 25, linetype = "dashed", color = "red") +
  geom_hline(yintercept = 50, linetype = "dashed", color = "blue") +

  scale_y_continuous(
    name = "Nilai RSE (%)",
    breaks = c(0, 25, 50),
    limits = c(-5, NA)
  ) +

  scale_color_manual(
    name = "Kategori RSE",
    limits = c("< 25%", "26–50%", "> 50%"),
    values = c(
      "< 25%"  = "#1B9E77",
      "26–50%" = "#FF8C00",
      "> 50%"  = "#7570B3"
    ),
    na.value = "#D73027",
    drop = FALSE
  ) +
  theme_minimal()

Satu kabupaten/kota memiliki nilai RSE = 0 sehingga diperlakukan sebagai data tidak tersedia (NA) dan ditampilkan dalam visualisasi sebagai titik berwarna merah.

write.xlsx(outputP,"BACKLOG1 2025_Prov.xlsx")
write.xlsx(output,"BACKLOG1 2025_Kako.xlsx")

5 Peta

library(sf)
petaku <- st_read("PetaSHP514_38.shp")

## Reading layer `PetaSHP514_38' from data source 
##   `D:\2. Pengembangan diri\1 Exercise Bagus\Housing Analysis\PetaSHP514_38.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 514 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.00971 ymin: -11.00766 xmax: 141.02 ymax: 6.076809
## Geodetic CRS:  WGS 84

petaku <- petaku %>%
  mutate(idkab = as.character(idkab))
hasil <- hasil %>%
  mutate(kabu = as.character(kabu))
names(petaku)

## [1] "fid"      "idkab"    "nmprov"   "nmkab"    "kdprov"   "kdkab"    "sumber"  
## [8] "periode"  "geometry"

kab_map <- petaku %>%
  left_join(
    hasil %>% select(kabu, RSE, RSE_kat),
    by = c("idkab" = "kabu")
  )
kab_map$RSE_kat <- factor(
  kab_map$RSE_kat,
  levels = c("< 25%", "26–50%", "> 50%")
)

---

Direktorat Statistik Kesejahteraan Rakyat, BPS, saptahas@bps.go.id