library(haven) # baca/tulis file Stata (.dta), fungsi labelled()
library(dplyr) # filter, mutate, group_by, summarise, joinHand-out: Konstruksi & Agregasi Data SUSENAS–PODES dengan R
Studi kasus: Balita Sakit (SUSENAS) & Posyandu (PODES) — data sintetis 119 kab/kota Pulau Jawa
ImportantDisklaimer: data sintetis, bukan data asli
Seluruh data pada hand-out ini adalah data sintetis — dibangkitkan komputer khusus untuk keperluan pelatihan. Struktur file, nama variabel (R101, R407, R1102, R705A, dst.), dan alur pengolahan meniru SUSENAS KOR Maret 2023 dan PODES 2021, tetapi seluruh angka, indikator, dan hasil regresi tidak menggambarkan kondisi riil kabupaten/kota mana pun dan bukan data resmi BPS. Jangan mengutip angka dari materi ini. Untuk data sesungguhnya, ajukan akses data mikro resmi BPS melalui silastik.bps.go.id.
1 Gambaran umum alur kerja
Analisis kabupaten/kota berbasis data mikro BPS selalu melewati tiga tahap:
- Konstruksi — dari file mentah (satu baris = satu individu/rumah tangga/desa) dibentuk dataset per variabel yang bersih: kode wilayah, variabel hasil recode, bobot, dan label.
- Agregasi — dataset per variabel diringkas menjadi satu angka per kabupaten/kota. Untuk survei (SUSENAS) wajib memakai bobot (
FWTuntuk individu,WERTuntuk rumah tangga); untuk sensus (PODES) cukup penjumlahan. - Penggabungan — semua indikator kabupaten di-merge dengan kunci
kode_kabmenjadi satu tabel basis analisis (119 baris untuk Pulau Jawa).
Hand-out ini mencontohkan satu variabel dari tiap sumber:
| Sumber | Variabel contoh | File mentah |
|---|---|---|
| SUSENAS | Balita punya masalah kesehatan (R1102) |
kor23_ind_1_sintetis.dta |
| PODES | Jumlah posyandu di desa (R705A) |
podes_desa_sintetis.dta |
2 Mengenal antarmuka RStudio
Sebelum mulai, kenali dulu empat panel utama RStudio — seluruh latihan akan bolak-balik di antara panel-panel ini:
Ringkasnya:
- Menu & Toolbar — membuka file, mengatur working directory (
Session > Set Working Directory), dan tombol Render untuk file Quarto (.qmd). - Source / Editor — tempat membuka dan mengedit skrip R (misalnya
11_stepbystep_quarto.R). Jalankan baris demi baris dengan Ctrl+Enter. - Console — hasil eksekusi kode tampil di sini; perintah juga bisa diketik langsung setelah tanda
>. - Environment — daftar objek yang sedang dimuat (data dan model); klik nama data untuk membuka tampilan tabelnya.
- Files / Plots / Packages / Help — grafik muncul di tab Plots; Files untuk navigasi folder proyek; Help untuk dokumentasi fungsi (mis. ketik
?lmdi Console).
Persiapan: dua paket ini yang dipakai sepanjang latihan.
3 SUSENAS — Konstruksi variabel “balita sakit”
3.1 Baca dan kenali data mentah
File mentah KOR individu: satu baris = satu anggota rumah tangga.
kor1 <- read_dta("Raw_Sintetis_DTA/kor23_ind_1_sintetis.dta")
dim(kor1)[1] 44702 12Kolom-kolom kunci untuk variabel ini:
head(kor1[, c("URUT", "R101", "R102", "R401", "R407", "R1102", "FWT")])| URUT | R101 | R102 | R401 | R407 | R1102 | FWT |
|---|---|---|---|---|---|---|
| 31010001 | 31 | 01 | 1 | 30 | 3 | 98.38 |
| 31010001 | 31 | 01 | 2 | 29 | 2 | 98.38 |
| 31010001 | 31 | 01 | 3 | 0 | 3 | 98.38 |
| 31010001 | 31 | 01 | 4 | 11 | 1 | 98.38 |
| 31010002 | 31 | 01 | 1 | 29 | 4 | 85.80 |
| 31010002 | 31 | 01 | 2 | 23 | 2 | 85.80 |
| Variabel | Arti |
|---|---|
URUT |
nomor urut rumah tangga (id ruta) |
R101 |
kode provinsi (31–36 di Jawa) |
R102 |
kode kabupaten/kota, 2 digit (character) |
R401 |
nomor urut anggota dalam ruta |
R407 |
umur (tahun) |
R1102 |
jenis keluhan kesehatan; kode 5 = punya masalah kesehatan |
FWT |
bobot individu (final weight) |
3.2 Bentuk kode kabupaten dan id individu
Kode BPS 4 digit dibentuk dengan menempelkan R101 dan R102:
kor1 <- kor1 %>%
mutate(
kode_kab = as.integer(paste0(R101, R102)), # mis. 33 + "74" -> 3374
kode_ind = paste0(URUT, R401) # id unik individu
)
kor1 %>% select(URUT, R101, R102, kode_kab, kode_ind) %>% head(4)| URUT | R101 | R102 | kode_kab | kode_ind |
|---|---|---|---|---|
| 31010001 | 31 | 01 | 3101 | 310100011 |
| 31010001 | 31 | 01 | 3101 | 310100012 |
| 31010001 | 31 | 01 | 3101 | 310100013 |
| 31010001 | 31 | 01 | 3101 | 310100014 |
Ini padanan egen kode_kab = concat(R101 R102) di Stata.
3.3 Filter balita dan recode variabel
Definisi: balita = umur di bawah 5 tahun (R407 < 5); sakit = keluhan kesehatan berkode 5 (R1102 == 5).
balita_sakit <- kor1 %>%
filter(R407 < 5) %>%
mutate(balita_sakit = as.integer(R1102 == 5)) %>%
select(URUT, FWT, kode_ind, kode_kab, balita_sakit)
nrow(balita_sakit)[1] 5293table(balita_sakit$balita_sakit)
0 1
3897 1396 3.4 Pasang label dan simpan dataset per variabel
Supaya file .dta tetap informatif saat dibuka di Stata/R, pasang value label dan variable label:
balita_sakit <- balita_sakit %>%
mutate(balita_sakit = labelled(balita_sakit,
labels = c(Tidak = 0L, Ya = 1L),
label = "Balita memiliki masalah kesehatan; 1=Ya, 0=Tidak"))
write_dta(balita_sakit, "Data_Sintetis_DTA/susenas_balita_sakit.dta")Tahap konstruksi selesai: satu file bersih, satu baris = satu balita.
3.5 Agregasi ke kabupaten dengan bobot FWT
SUSENAS adalah survei — setiap responden “mewakili” sejumlah orang di populasi sebesar bobotnya. Proporsi kabupaten dihitung dengan weighted.mean():
kab_balita <- balita_sakit %>%
group_by(kode_kab) %>%
summarise(
n_balita = n(), # jumlah sampel, tanpa bobot
pct_balita_sakit = weighted.mean(as.numeric(balita_sakit),
w = FWT) * 100,
.groups = "drop")
head(kab_balita)| kode_kab | n_balita | pct_balita_sakit |
|---|---|---|
| 3101 | 34 | 20.218433 |
| 3171 | 40 | 8.874837 |
| 3172 | 45 | 10.748338 |
| 3173 | 54 | 14.347336 |
| 3174 | 28 | 20.517733 |
| 3175 | 33 | 25.674686 |
Mengapa bobot wajib? Bandingkan hasil tanpa vs dengan bobot:
balita_sakit %>%
group_by(kode_kab) %>%
summarise(
tanpa_bobot = mean(as.numeric(balita_sakit)) * 100,
dengan_bobot = weighted.mean(as.numeric(balita_sakit), w = FWT) * 100
) %>% head(5)| kode_kab | tanpa_bobot | dengan_bobot |
|---|---|---|
| 3101 | 20.58824 | 20.218433 |
| 3171 | 10.00000 | 8.874837 |
| 3172 | 11.11111 | 10.748338 |
| 3173 | 18.51852 | 14.347336 |
| 3174 | 25.00000 | 20.517733 |
Selisihnya bisa beberapa poin persen — desain sampling SUSENAS tidak self-weighting, sehingga rata-rata tanpa bobot bias.
Tip
Untuk standard error yang benar (bukan sekadar titik estimasi), gunakan paket survey: svydesign(id = ~1, weights = ~FWT, data = ...) lalu svyby(~balita_sakit, ~kode_kab, design, svymean).
4 PODES — Akumulasi fasilitas kesehatan
4.1 Baca dan kenali data mentah desa
PODES adalah sensus seluruh desa/kelurahan — satu baris = satu desa, tanpa bobot.
podes <- read_dta("Raw_Sintetis_DTA/podes_desa_sintetis.dta")
dim(podes)[1] 3051 18head(podes[, c("R101", "R102", "R103", "R104", "R705A")])| R101 | R102 | R103 | R104 | R705A |
|---|---|---|---|---|
| 31 | 01 | 001 | 001 | 2 |
| 31 | 01 | 001 | 002 | 5 |
| 31 | 01 | 001 | 003 | 3 |
| 31 | 01 | 001 | 004 | 3 |
| 31 | 01 | 001 | 005 | 1 |
| 31 | 01 | 001 | 006 | 6 |
| Variabel | Arti |
|---|---|
R101–R104 |
kode provinsi, kab/kota, kecamatan, desa |
R704*K2 |
jumlah tiap jenis fasilitas kesehatan (blok 704) |
R705A |
jumlah posyandu di desa |
4.2 Akumulasi ke kabupaten
Karena sensus, agregasi cukup sum() semua desa dalam kabupaten:
kab_posyandu <- podes %>%
mutate(kode_kab = as.integer(paste0(R101, R102))) %>%
group_by(kode_kab) %>%
summarise(
n_desa = n(),
jumlah_posyandu = sum(R705A, na.rm = TRUE),
.groups = "drop")
head(kab_posyandu)| kode_kab | n_desa | jumlah_posyandu |
|---|---|---|
| 3101 | 25 | 76 |
| 3171 | 33 | 110 |
| 3172 | 26 | 112 |
| 3173 | 30 | 118 |
| 3174 | 15 | 57 |
| 3175 | 23 | 81 |
4.3 Jadikan rasio agar sebanding antar kabupaten
Jumlah absolut menguntungkan kabupaten besar. Untuk perbandingan yang adil, bagi dengan populasi (di sini populasi sintetis dari tabel referensi):
ref <- read_dta("Data_Sintetis_DTA/ref_kabupaten_jawa.dta")
kab_posyandu <- kab_posyandu %>%
left_join(ref %>% select(kode_kab, nama_kab, populasi), by = "kode_kab") %>%
mutate(posyandu_per_100k = jumlah_posyandu / as.numeric(populasi) * 1e5)
kab_posyandu %>%
select(kode_kab, nama_kab, n_desa, jumlah_posyandu, posyandu_per_100k) %>%
head()| kode_kab | nama_kab | n_desa | jumlah_posyandu | posyandu_per_100k |
|---|---|---|---|---|
| 3101 | Kepulauan Seribu | 25 | 76 | 2.364424 |
| 3171 | Jakarta Selatan | 33 | 110 | 9.173164 |
| 3172 | Jakarta Timur | 26 | 112 | 11.364397 |
| 3173 | Jakarta Pusat | 30 | 118 | 5.479991 |
| 3174 | Jakarta Barat | 15 | 57 | 4.107276 |
| 3175 | Jakarta Utara | 23 | 81 | 12.777013 |
5 Gabungkan menjadi basis data regresi
Kunci penggabungan selalu kode_kab:
data_kab <- kab_balita %>%
left_join(kab_posyandu, by = "kode_kab") %>%
select(kode_kab, nama_kab, pct_balita_sakit,
jumlah_posyandu, posyandu_per_100k)
nrow(data_kab)[1] 119head(data_kab)| kode_kab | nama_kab | pct_balita_sakit | jumlah_posyandu | posyandu_per_100k |
|---|---|---|---|---|
| 3101 | Kepulauan Seribu | 20.218433 | 76 | 2.364424 |
| 3171 | Jakarta Selatan | 8.874837 | 110 | 9.173164 |
| 3172 | Jakarta Timur | 10.748338 | 112 | 11.364397 |
| 3173 | Jakarta Pusat | 14.347336 | 118 | 5.479991 |
| 3174 | Jakarta Barat | 20.517733 | 57 | 4.107276 |
| 3175 | Jakarta Utara | 25.674686 | 81 | 12.777013 |
Visual cepat hubungan kedua variabel:
plot(data_kab$posyandu_per_100k, data_kab$pct_balita_sakit,
xlab = "Posyandu per 100.000 penduduk",
ylab = "% balita sakit (tertimbang FWT)",
pch = 19, col = "steelblue")
abline(lm(pct_balita_sakit ~ posyandu_per_100k, data = data_kab),
col = "firebrick", lwd = 2)
Dan regresi sederhananya:
summary(lm(pct_balita_sakit ~ posyandu_per_100k, data = data_kab))
Call:
lm(formula = pct_balita_sakit ~ posyandu_per_100k, data = data_kab)
Residuals:
Min 1Q Median 3Q Max
-17.3897 -6.4269 -0.4881 4.9274 27.0829
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.6510 1.1375 25.188 < 2e-16 ***
posyandu_per_100k -0.3923 0.1365 -2.873 0.00483 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.179 on 117 degrees of freedom
Multiple R-squared: 0.06589, Adjusted R-squared: 0.05791
F-statistic: 8.253 on 1 and 117 DF, p-value: 0.0048316 Analisa OLS dengan variabel lengkap
Bagian sebelumnya hanya memakai dua variabel contoh. Dengan pola yang sama, seluruh variabel lain dapat dikonstruksi (skrip 03–09) lalu digabungkan (skrip 10) menjadi basis data Hasil_Agregasi_DTA/data_kabupaten_jawa.dta — 119 kabupaten/kota dengan indikator SUSENAS dan PODES lengkap. Bagian ini memakai basis data tersebut untuk analisa OLS.
6.1 Baca basis data dan siapkan variabel
data_jawa <- read_dta("Hasil_Agregasi_DTA/data_kabupaten_jawa.dta") %>%
mutate(kota = as.integer(tipe == 2)) # dummy: 1 = Kota, 0 = Kabupaten
dim(data_jawa)[1] 119 16names(data_jawa) [1] "kode_kab" "nama_kab" "provinsi"
[4] "tipe" "populasi" "pct_balita_sakit"
[7] "pct_imunisasi" "pct_keluarga_merokok" "pct_ortu_merokok"
[10] "mean_pendidikan_ibu" "mean_pengeluaran_susu" "jumlah_posyandu"
[13] "posyandu_per_100k" "total_faskes" "faskes_per_100k"
[16] "kota" Model yang disarankan — determinan kesehatan balita antar kabupaten:
\[\text{pct\_balita\_sakit}_i = \beta_0 + \beta_1\,\text{pct\_imunisasi}_i + \beta_2\,\text{pct\_keluarga\_merokok}_i + \beta_3\,\text{mean\_pendidikan\_ibu}_i + \beta_4\,\text{faskes\_per\_100k}_i + \beta_5\,\text{kota}_i + \varepsilon_i\]
Ekspektasi tanda koefisien: imunisasi (−), keluarga merokok (+), pendidikan ibu (−), faskes per kapita (−).
6.2 Eksplorasi sebelum regresi
Selalu mulai dari statistik deskriptif dan korelasi:
data_jawa %>%
select(pct_balita_sakit, pct_imunisasi, pct_keluarga_merokok,
mean_pendidikan_ibu, faskes_per_100k, kota) %>%
summary() pct_balita_sakit pct_imunisasi pct_keluarga_merokok mean_pendidikan_ibu
Min. : 8.875 Min. :29.90 Min. :32.36 Min. :1.434
1st Qu.:20.321 1st Qu.:48.03 1st Qu.:48.53 1st Qu.:2.283
Median :25.655 Median :54.91 Median :53.91 Median :2.618
Mean :26.452 Mean :55.98 Mean :53.68 Mean :2.601
3rd Qu.:32.222 3rd Qu.:64.33 3rd Qu.:59.98 3rd Qu.:2.876
Max. :55.712 Max. :82.76 Max. :72.83 Max. :4.568
faskes_per_100k kota
Min. : 1.155 Min. :0.0000
1st Qu.: 3.948 1st Qu.:0.0000
Median : 6.374 Median :0.0000
Mean : 9.282 Mean :0.2857
3rd Qu.: 9.525 3rd Qu.:1.0000
Max. :67.993 Max. :1.0000 Cara membaca: persentase balita sakit sangat bervariasi antar kabupaten — dari sekitar 9% sampai 56%, dengan rata-rata ±26%. Variasi yang lebar ini justru bagus untuk regresi: ada “sesuatu” yang perlu dijelaskan. Cakupan imunisasi juga beragam (30–83%), begitu pula rasio faskes (1–68 per 100 ribu penduduk — nilai maksimum yang jauh dari kuartil 3 menandakan ada outlier, kemungkinan kota kecil berpenduduk sedikit).
round(cor(data_jawa %>%
select(pct_balita_sakit, pct_imunisasi, pct_keluarga_merokok,
mean_pendidikan_ibu, faskes_per_100k)), 2) pct_balita_sakit pct_imunisasi pct_keluarga_merokok
pct_balita_sakit 1.00 -0.45 0.24
pct_imunisasi -0.45 1.00 -0.15
pct_keluarga_merokok 0.24 -0.15 1.00
mean_pendidikan_ibu -0.52 0.61 -0.16
faskes_per_100k -0.24 0.39 -0.14
mean_pendidikan_ibu faskes_per_100k
pct_balita_sakit -0.52 -0.24
pct_imunisasi 0.61 0.39
pct_keluarga_merokok -0.16 -0.14
mean_pendidikan_ibu 1.00 0.23
faskes_per_100k 0.23 1.00Cara membaca: baris/kolom pertama menunjukkan korelasi outcome dengan tiap prediktor — semuanya searah dugaan: negatif dengan imunisasi (−0,45) dan pendidikan ibu (−0,52), positif dengan keluarga merokok (+0,24). Perhatikan juga korelasi antar prediktor: imunisasi × pendidikan ibu mencapai 0,61 — sinyal awal multikolinearitas yang nanti kita konfirmasi dengan VIF.
6.3 Visualisasi deskriptif
Angka ringkasan perlu didampingi gambar — sebaran, perbandingan antar kelompok, dan hubungan antar variabel jauh lebih cepat terbaca secara visual.
(a) Sebaran outcome. Histogram memperlihatkan bentuk sebaran % balita sakit antar 119 kabupaten/kota (garis putus-putus = rata-rata):
hist(data_jawa$pct_balita_sakit,
breaks = 20, col = "steelblue", border = "white",
main = "Sebaran % Balita Sakit antar Kabupaten/Kota",
xlab = "% balita sakit (tertimbang FWT)", ylab = "Frekuensi")
abline(v = mean(data_jawa$pct_balita_sakit), col = "firebrick",
lwd = 2, lty = 2)
Cara membaca: sebaran cenderung simetris dengan puncak di sekitar 20–30% dan ekor kanan yang memanjang — beberapa kabupaten punya persentase balita sakit jauh di atas rata-rata (garis putus-putus merah). Tidak ada indikasi sebaran ekstrem yang menuntut transformasi (mis. log) pada outcome.
(b) Perbandingan antar provinsi. Boxplot menunjukkan variasi antar dan di dalam provinsi:
par(mar = c(9, 4, 3, 1))
boxplot(pct_balita_sakit ~ provinsi, data = data_jawa,
col = "lightsteelblue", las = 2, xlab = "",
main = "% Balita Sakit menurut Provinsi",
ylab = "% balita sakit")
par(mar = c(5, 4, 4, 2))Cara membaca: median antar provinsi berbeda nyata — DKI Jakarta terendah (median ±17%) dan Banten tertinggi (±34%), selisih hampir dua kali lipat. Kotak yang lebar (mis. Jawa Barat) menandakan keragaman antar kabupaten di dalam provinsi itu juga besar. Perbedaan level antar provinsi seperti ini yang membuat dummy wilayah sering dimasukkan ke model.
(c) Kabupaten vs kota. Apakah wilayah perkotaan berbeda?
boxplot(pct_balita_sakit ~ kota, data = data_jawa,
names = c("Kabupaten", "Kota"),
col = c("lightsteelblue", "lightsalmon"),
main = "% Balita Sakit: Kabupaten vs Kota",
ylab = "% balita sakit")
Cara membaca: kota sedikit lebih rendah daripada kabupaten (rata-rata ±24% vs ±27%), tetapi kedua kotak masih banyak tumpang tindih — perbedaannya tidak dramatis. Ini konsisten dengan hasil regresi nanti: koefisien dummy kota bertanda negatif tetapi tidak signifikan.
(d) Matriks scatterplot. Semua pasangan variabel sekaligus — baris pertama memperlihatkan hubungan outcome dengan tiap prediktor, sel-sel lainnya memperlihatkan korelasi antar prediktor (bahan diskusi multikolinearitas):
pairs(data_jawa[, c("pct_balita_sakit", "pct_imunisasi",
"pct_keluarga_merokok", "mean_pendidikan_ibu",
"faskes_per_100k")],
labels = c("% balita\nsakit", "% imunisasi", "% keluarga\nmerokok",
"pendidikan\nibu", "faskes\nper 100k"),
pch = 19, col = adjustcolor("steelblue", 0.6),
main = "Matriks Scatterplot: Outcome vs Prediktor")
Cara membaca: fokus ke baris pertama — pola menurun terlihat pada imunisasi dan pendidikan ibu (hubungan negatif dengan outcome), pola menaik samar pada keluarga merokok. Di luar baris pertama, panel imunisasi × pendidikan ibu memperlihatkan pola menaik yang cukup rapat: inilah wujud visual korelasi 0,61 antar prediktor tadi.
6.4 Model 1: regresi sederhana
ols1 <- lm(pct_balita_sakit ~ pct_imunisasi, data = data_jawa)
summary(ols1)
Call:
lm(formula = pct_balita_sakit ~ pct_imunisasi, data = data_jawa)
Residuals:
Min 1Q Median 3Q Max
-18.1976 -5.2648 -0.3567 4.3136 25.4558
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 47.37559 3.87524 12.225 < 2e-16 ***
pct_imunisasi -0.37376 0.06782 -5.511 2.16e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.462 on 117 degrees of freedom
Multiple R-squared: 0.2061, Adjusted R-squared: 0.1993
F-statistic: 30.37 on 1 and 117 DF, p-value: 2.157e-07Cara membaca: setiap kenaikan 1 poin persen cakupan imunisasi berasosiasi dengan penurunan % balita sakit sebesar ±0,37 poin (p < 0,001). R² = 0,21 artinya imunisasi sendirian “menjelaskan” sekitar 21% variasi antar kabupaten. Tetapi ingat: model satu variabel rawan omitted variable bias — koefisien ini masih “menampung” efek variabel lain yang berkorelasi dengan imunisasi (terutama pendidikan ibu).
6.5 Model 2: regresi berganda
ols2 <- lm(pct_balita_sakit ~ pct_imunisasi + pct_keluarga_merokok +
mean_pendidikan_ibu + faskes_per_100k + kota,
data = data_jawa)
summary(ols2)
Call:
lm(formula = pct_balita_sakit ~ pct_imunisasi + pct_keluarga_merokok +
mean_pendidikan_ibu + faskes_per_100k + kota, data = data_jawa)
Residuals:
Min 1Q Median 3Q Max
-22.1456 -4.8483 -0.3046 4.8743 25.5115
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.34234 6.73183 6.587 1.48e-09 ***
pct_imunisasi -0.16284 0.08533 -1.908 0.058882 .
pct_keluarga_merokok 0.16626 0.08666 1.918 0.057576 .
mean_pendidikan_ibu -6.43909 1.81798 -3.542 0.000578 ***
faskes_per_100k -0.03965 0.08229 -0.482 0.630839
kota -2.04015 1.67600 -1.217 0.226036
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.895 on 113 degrees of freedom
Multiple R-squared: 0.3325, Adjusted R-squared: 0.303
F-statistic: 11.26 on 5 and 113 DF, p-value: 8.012e-09confint(ols2) 2.5 % 97.5 %
(Intercept) 31.005363102 57.67931894
pct_imunisasi -0.331885358 0.00621379
pct_keluarga_merokok -0.005435889 0.33796023
mean_pendidikan_ibu -10.040830900 -2.83734238
faskes_per_100k -0.202676929 0.12337538
kota -5.360613229 1.28030744Cara membaca: semua tanda koefisien sesuai ekspektasi. Prediktor terkuat adalah pendidikan ibu: naik satu jenjang rata-rata pendidikan ibu berasosiasi dengan turunnya % balita sakit ±6,4 poin (p < 0,001) — konsisten dengan temuan umum literatur kesehatan anak. Imunisasi (−0,16) dan keluarga merokok (+0,17) berada di ambang signifikansi (p ≈ 0,06), sedangkan faskes per kapita dan dummy kota tidak signifikan setelah variabel lain dikendalikan. Interval kepercayaan (confint) yang memuat nol menandai koefisien yang tidak signifikan pada taraf 5%.
NoteMomen didaktik: omitted variable bias
Koefisien pct_imunisasi menyusut dari -0.37 (Model 1) menjadi -0.16 (Model 2). Sebabnya: imunisasi berkorelasi kuat dengan pendidikan ibu, sehingga pada model sederhana sebagian efek pendidikan “menumpang” pada koefisien imunisasi. Inilah ilustrasi klasik omitted variable bias.
6.6 Diagnostik OLS
Empat plot diagnostik standar:
par(mfrow = c(2, 2))
plot(ols2)
par(mfrow = c(1, 1))Cara membaca: Residuals vs Fitted tidak menunjukkan pola melengkung (hubungan linear memadai); Q-Q plot mengikuti garis diagonal dengan sedikit simpangan di ekor (normalitas residual dapat diterima untuk n = 119); Scale-Location relatif datar (varians residual stabil); dan di Residuals vs Leverage tidak ada titik yang melewati jarak Cook — tidak ada observasi tunggal yang “menyetir” hasil regresi.
Multikolinearitas (VIF) — paket car (pasang dengan install.packages("car") bila belum ada):
library(car)
vif(ols2) pct_imunisasi pct_keluarga_merokok mean_pendidikan_ibu
1.818361 1.046816 1.651193
faskes_per_100k kota
1.263748 1.094486 Cara membaca: semua VIF di bawah 2 (tertinggi imunisasi ±1,8) — jauh dari ambang masalah (rule of thumb: VIF > 10 serius, > 5 perlu waspada). Jadi walau imunisasi dan pendidikan ibu berkorelasi 0,61, kolinearitasnya belum mengganggu presisi estimasi; kedua variabel aman dipertahankan dalam satu model.
Heteroskedastisitas (uji Breusch-Pagan) — paket lmtest:
library(lmtest)
bptest(ols2)
studentized Breusch-Pagan test
data: ols2
BP = 7.2642, df = 5, p-value = 0.2017Cara membaca: p-value = 0,20 > 0,05 → gagal tolak H0 varians residual konstan (homoskedastis). Standard error OLS biasa dapat dipakai; robust standard error (paket sandwich, vcovHC) tidak wajib di sini — tetapi tetap praktik yang baik untuk data riil.
6.7 Perbandingan model
cat("R-squared model 1 :", round(summary(ols1)$r.squared, 3), "\n")R-squared model 1 : 0.206 cat("R-squared model 2 :", round(summary(ols2)$r.squared, 3), "\n")R-squared model 2 : 0.333 Cara membaca: R² naik dari 0,21 menjadi 0,33 — penambahan empat variabel menjelaskan tambahan ±13 poin persen variasi outcome. Sisanya (±67%) tak terjelaskan model: wajar untuk data lintas-wilayah, dan menjadi pengingat bahwa R² tinggi bukan tujuan utama — yang penting spesifikasi masuk akal dan koefisien dapat diinterpretasikan.
Warning
Ingat: data ini sintetis — tanda dan besaran koefisien mengikuti desain generator data, bukan temuan empiris. Gunakan hanya untuk memahami alur kerja dan interpretasi output.
7 Rangkuman dan latihan
| Tahap | SUSENAS (survei) | PODES (sensus desa) |
|---|---|---|
| Unit data mentah | individu / rumah tangga | desa |
| Kode wilayah | paste0(R101, R102) |
paste0(R101, R102) |
| Konstruksi | filter() + recode + labelled() |
pilih kolom faskes |
| Agregasi | weighted.mean(x, w = FWT) |
sum(x) |
| Hasil per kabupaten | persentase / rata-rata | jumlah / rasio per penduduk |
Latihan mandiri — ulangi seluruh alur untuk variabel lain:
- SUSENAS: imunisasi Hepatitis B (
R1404A_TGLterisi, filekor23_ind_2_sintetis.dta) → % balita terimunisasi per kabupaten. - SUSENAS: keluarga merokok (
R1206/R1207) — perhatikan agregasi maksimum per rumah tangga sebelum dihitung proporsinya. - PODES: jumlah puskesmas (
R704CK2 + R704DK2) per kabupaten, lalu jadikan rasio per 100.000 penduduk. - OLS: ganti outcome menjadi
pct_imunisasi— prediktor apa yang masuk akal? Coba juga tambahkanlog(mean_pengeluaran_susu)sebagai prediktor dan bandingkan hasilnya.
Skrip lengkap seluruh variabel tersedia di folder Script_R/ (00_ref_kabupaten_jawa.R sampai 11_stepbystep_quarto.R — skrip 11 adalah versi step-by-step yang persis mengikuti hand-out ini).
Daftar Pustaka
- Badan Pusat Statistik. (2023). Pedoman Pencacahan Survei Sosial Ekonomi Nasional (SUSENAS) Maret 2023. Jakarta: BPS.
- Badan Pusat Statistik. (2021). Pedoman Pendataan Potensi Desa (PODES) 2021. Jakarta: BPS.
- Fox, J., & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). Thousand Oaks, CA: Sage. (paket
car) - Lumley, T. (2010). Complex Surveys: A Guide to Analysis Using R. Hoboken, NJ: John Wiley & Sons. (paket
survey) - R Core Team. (2024). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. https://www.R-project.org/
- Wickham, H., François, R., Henry, L., Müller, K., & Vaughan, D. (2023). dplyr: A Grammar of Data Manipulation. R package. https://dplyr.tidyverse.org
- Wickham, H., & Miller, E. (2023). haven: Import and Export ‘SPSS’, ‘Stata’ and ‘SAS’ Files. R package. https://haven.tidyverse.org
- Wooldridge, J. M. (2020). Introductory Econometrics: A Modern Approach (7th ed.). Boston, MA: Cengage Learning.
- Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7–10. (paket
lmtest)
Pengingat terakhir: seluruh materi ini menggunakan data sintetis untuk pelatihan — bukan data resmi BPS dan tidak menggambarkan kondisi riil wilayah mana pun.