Load Packages

pacman::p_load(foreign, haven, car, dplyr, tidyr, survey, janitor, magrittr, ggplot2)

Setting Working Directory

setwd("D:\\Materi\\")

Load Dataset

kor23gab <- read_sav("KOR23GAB_2609.sav")

Selected Unit Analysis

Exclusive Breastfeeding (EBF) atau ASI Eksklusif ditujukan untuk bayi dengan usia dibawah 6 bulan. Oleh karena itu dalam penghitungan perlu dilakukan penyeleksian unit analisis terlebih dahulu.

P05 <- filter(kor23gab, R1401<=5)

Manipulating working variables

Untuk meringankan beban kerja dari proses penghitungan, perlu dilakukan filter variabel yang digunakan. Terkait indikator EBF, variabel yang digunakan antara lain: kode provinsi, kode kab/kota, kode renumbering NKS, kode renumbering no urut ruta, strata, klasifikasi kota desa, umur (bulan), variabel pemberian makan baduta udan penimbang.

dataku <- P05 %>%
  select(R101, R102, WI1, WI2, R105, R1401, R1417A:R1420H, R1421, FWT)
View(dataku)

Calculate the indicator

Exclusive Breastfeeding (EBF) atau ASI Eksklusif adalah air susu ibu yang diberikan kepada bayi sejak dilahirkan selama 6 (enam) bulan (0–5 bulan 29 hari), tanpa menambahkan dan/atau mengganti dengan makanan atau minuman lain kecuali obat, vitamin, dan mineral dalam 24 jam terakhir.

WHO Global Strategy for Infant and Young Child Feeding merekomendasikan agar bayi disusui secara eksklusif sampai mereka berusia enam bulan. ASI Eksklusif adalah pilihan teraman dan tersehat untuk anak-anak, menjamin bayi mendapatkan sumber makanan yang aman, bersih, sehat, dan mudah diakses.

  1. Merupakan indikator kinerja utama pada sektor K/L

  2. Sebagai bagian dari intervensi spesifik sekaligus faktor determinan penanganan stunting

library(dplyr)
dataku <- dataku %>%
  mutate(kabu = R101*100 + R102) %>%
  rename(prov="R101", kab="R102") %>%
  mutate(psu=prov*10000000 + kab*100000 + WI1, 
         denom=100, 
         strata= prov*10000000 + kab*100000 + R105,
         theta = dplyr::case_when(dataku$R1417A==1 & dataku$R1418==1 & 
                                  dataku$R1419A==5 & dataku$R1419B==5 & dataku$R1419C==5 &
                                  dataku$R1420A==5 & dataku$R1420B==5 & dataku$R1420C==5 &
                                  dataku$R1420D==5 & dataku$R1420E==5 & dataku$R1420F==5 &
                                  dataku$R1420G==5 & dataku$R1420H==5 & dataku$R1421==5~100,)) %>%
  mutate(theta = ifelse(is.na(theta), 0, theta))

Set Design Sampling

susenas.design<- svydesign(id=~psu, strata=~strata, 
                            data = dataku, weights=~FWT)
summary(susenas.design$prob)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0002617 0.0033376 0.0069436 0.0197024 0.0161337 0.4777765

Set Output Statistics

options(survey.adjust.domain.lonely=TRUE)
options(survey.lonely.psu="adjust")
hasil = svyby(formula = ~theta, 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$`theta/denom`*100,2)
hasil$SE = round(hasil$`se.theta/denom`*100,2)
hasil$VAR = round(hasil$SE*hasil$SE,2)
hasil$CI_LOWER = round(hasil$`ci_l`*100,2)
hasil$CI_UPPER = round(hasil$`ci_u`*100,2)
hasil$RSE = round(hasil$`cv%`,2)
hasil$DEFF = round(hasil$DEff,2)
View(hasil)

District Level Estimate

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")
(output[order(output$Kako, decreasing = FALSE), ]   )

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\le50%\), estimasi perlu dilakukan dengan hati-hati

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

ggplot(data = hasil, aes(x=hasil$kabu, y=hasil$RSE)) +
  geom_line(colour = "orange") +
  labs(title = "Kabupaten/Kota vs Nilai RSE")

Dari plot di atas, tampak nilai RSE>25% hampir ditemukan di setiap wilayah. Oleh karena itu perlu dilakukan pengestimasian melalui cara lain, misalnya melalui Small Area Estimation (SAE).

Export to CSV

write.csv(output,"RSE_DIRECT_EBF_KAKO_2023.csv")

rm(list=ls()) #Clear environment

References

  1. BPS. 2023. Profil Statistik Kesehatan 2023. Jakarta: Badan Pusat Statistik.
  2. Confidentiality and relative standard error, Australian Bureau Statistics. https://www.abs.gov.au/statistics/microdata-tablebuilder/tablebuilder/confidentiality-and-relative-standard-error
  3. Indicators for assessing infant and young child feeding practices: definitions and measurement methods. 2021. Geneva: World Health Organization and the United Nations Children’s Fund (UNICEF).

Direktorat Statistik Kesejahteraan Rakyat, BPS,