projek3

Masukkan Password untuk Melanjutkan

Submit

0.1 library

library(dplyr)
library(ggplot2)
library(lmtest)
library(lme4) 
library(car)
library(tseries)
library(ResourceSelection)
library(psych)
library(foreign)
library(WRS2)
library(pROC)
library(nortest)
library(survey)

1 data preparation

1.1 data provinsi sumatera barat

Menyimpan data susenas hanya untuk wilayah sumatera barat

https://docs.google.com/document/d/1_JLASE7oz14lXtFWylf1Vl05Wf-4nameoMuxNyIGbvY/edit?usp=sharing

df1a <- read.dbf("~/Downloads/Survei-Sosial-Ekonomi-Nasional-2024-Maret-KOR/811 Susenas 202403 Kor/ssn202403_kor_ind1.dbf")
df2a <- read.dbf("~/Downloads/Survei-Sosial-Ekonomi-Nasional-2024-Maret-KOR/811 Susenas 202403 Kor/ssn202403_kor_ind2.dbf")
df3a <- read.dbf("~/Downloads/Survei-Sosial-Ekonomi-Nasional-2024-Maret-KOR/811 Susenas 202403 Kor/ssn202403_kor_rt.dbf")
df4a <- read.dbf("~/Downloads/Survei-Sosial-Ekonomi-Nasional-2024-Maret-Modul-Konsumsi-dan-Pengeluaran/812 Susenas 202403 KP/ssn202403_kp_blok43.dbf")

# tiap data frame simpan untuk provinsi sumatera barat Saja
df1 <- df1a %>% filter(R101==13)  #individu
df2 <- df2a %>% filter(R101==13)  #individu
df3 <- df3a %>% filter(R101==13)  #ruta
df4 <- df4a %>% filter(R101==13)  #ruta kemiskinan

1.2 data aggregate susenas KOR

menggabungkan data Susenas KOR berdasarkan suatu kunci dan unit analisisnya
hanya Kepala Rumah Tangga (KRT)

# Cari nama kolom yang sama di ketiga data frame
common_id <- Reduce(intersect, list(names(df1), names(df2)))

data <- inner_join(df1,df2, by = common_id)

## untuk KRT
data <- data %>% filter(R403==1)

common_id2 <- Reduce(intersect, list(names(data), names(df3)))
data_f <- data %>% 
  inner_join(df3, by = common_id2)

## mengecek bahwa data nya match
nrow(data_f) == nrow(data)  # Harus TRUE
nrow(data_f)

1.3 data aggregat susenas modul dengan GK

menggabungkan data susenas kor dengan modul pengeluaran termasuk garis kemiskinan kabupaten/kota untuk mengetahui penduduk yang miskin/tidak

# Membuat data frame yaitu kabupaten dan kota beserta garis kemiskinannya
KABU_GK <- data.frame(
  R102 = c(1:12, 71:77),
  GK = c(461088, 589787, 569605, 565067, 554665, 580392, 
         558963, 582794, 515608, 551340, 622190, 657354,
         736786, 569869, 524196, 634303, 658640, 648230, 609286)
)

data_k <- inner_join(KABU_GK,df4, by = "R102")

# Cari nama kolom yang sama di ketiga data frame
common_id3 <- Reduce(intersect, list(names(data_f), names(data_k)))
common_id3

data_f <- data_f %>% 
  inner_join(data_k, by = common_id3)

## mengecek bahwa data nya match
nrow(data_k) == nrow(data_f)
nrow(data_f)

1.4 unit dan variabel analisis dan melihat NA

unit analisis adalah bekerja/sementara tidak bekerja

# untuk memastikan bahwa NA di R612 (jenjang pendidikan tertinggi) karena   
# R610 adalah tidak bersekolah, lalu yg tidak bersekolah di recode 0 (<SMA)
sum(is.na(data_f$R612))
sum(data_f$R610==1)
all(is.na(data_f$R612[data_f$R610 == 1]))
data_f$R612[is.na(data_f$R612)] <- 0

# data hanya memuat var yang diteliti dan cek apakah ada NA
data_f %>% select(
  R105, R1802, R301, R405, R407, R612,R706, R703_A,R704,
  R2207, R2209A,R2209B, R2209C, R708, R1102,
  R1701, R1702, R1703, R1704, R1705, R1706, R1707, R1708,
  R102,R705,PSU, SSU, FWT, URUT, WI1, WI2,GK,KAPITA) %>%
  filter(R703_A == "A" | (R703_A != "A" & R705 == 1)) %>%
  summarise(across(everything(), ~sum(is.na(.))))

data_f1 <- data_f %>% select(
  R105, R1802, R301, R405, R407, R612,R706, R703_A,R704,
  R2207, R2209A,R2209B, R2209C, R708, R1102,
  R1701, R1702, R1703, R1704, R1705, R1706, R1707, R1708,
  R102,R705,PSU, SSU, FWT, URUT, WI1, WI2,GK,KAPITA,STRATA) %>%
  filter(R703_A == "A" | (R703_A != "A" & R705 == 1))

nrow(data_f1)

1.5 recode variabel yang bersesuain

termausk memberikan nama pada variabel yang dianalisis dan memberikan kategori terentu.

# Pilih variabel-variabel kategorik yang dibutuhkan
data1 <- data_f1 %>%
  mutate(
    # 1. Variabel Karakteristik
    perdesaan = factor(
      if_else(R105 == 2, 1L, 0L),
      levels = c(0, 1),
      labels = c("Perkotaan", "Perdesaan")
    ),
    status_rumah = factor(
      if_else(R1802 %in% c(1,4), 0L, 1L),
      levels = c(0, 1),
      labels = c("Milik Sendiri/Dinas", "Selainnya")
    ),
    art_banyak = factor(
      if_else(R301 > 4, 1L, 0L),
      levels = c(0, 1),
      labels = c("ART ≤4", "ART >4")
    ),
    perempuan = factor(
      if_else(R405 == 2, 1L, 0L),
      levels = c(0, 1),
      labels = c("Laki-laki", "Perempuan")
    ),
    usia_tua = factor(
      if_else(R407 > 45, 1L, 0L),
      levels = c(0, 1),
      labels = c("Usia ≤45", "Usia >45")
    ),
    pendidikan_tinggi = factor(
      if_else(R612 > 17, 1L, 0L),
      levels = c(0, 1),
      labels = c("≤SMA", ">SMA")
    ),
    pertanian = factor(
      if_else(R706 %in% 1:6, 1L, 0L),
      levels = c(0, 1),
      labels = c("NonPertanian", "Pertanian")
    ),
    miskin = factor(
      if_else(KAPITA < GK , 1L, 0L),
      levels = c(0, 1),
      labels = c("Tidak miskin", "Miskin")
    ),
    
     # 2. Variabel Kebijakan
    dapat_bpnt = factor(
      if_else(R2207 == 1, 1L, 0L),
      levels = c(0, 1),
      labels = c("Tidak Dapat BPNT", "Dapat BPNT")
    ),
    dapat_blt = factor(
      if_else(R2209A == 1, 1L, 0L),
      levels = c(0, 1),
      labels = c("Tidak Dapat BLT", "Dapat BLT")
    ),
    dapat_pktd = factor(
      if_else(R2209B == 1, 1L, 0L),
      levels = c(0, 1),
      labels = c("Tidak Dapat PKTD", "Dapat PKTD")
    ),
    dapat_bp = factor(
      if_else(R2209C == 1, 1L, 0L),
      levels = c(0, 1),
      labels = c("Tidak Dapat Bantuan Pangan", "Dapat Bantuan Pangan")
    ),
    # tingkat intensitas pekerjaan
    keluhan = factor(
      if_else(R1102 == 1, 1L, 0L),
      levels = c(0, 1),
      labels = c("Tidak ada", "Ada")),
    jam_kerja = R708, #jumlah jam kerja utama
      
    # Variabel FIES
    fies1 = R1701,          # Khawatir tidak cukup makanan
    fies2 = R1702,          # Tidak makan makanan bergizi
    fies3 = R1703,          # Makan sedikit jenis makanan
    fies4 = R1704,          # Melewatkan satu makan
    fies5 = R1705,          # Makan lebih sedikit 
    fies6 = R1706,          # Kehabisan makanan
    fies7 = R1707,          # Lapar tapi tidak makan
    fies8 = R1708,          # Tidak makan selama sehari penuh
    
    kab_kota = R102        # kabupaten/kota
    
  ) %>% select(
    # Sampling design variables
    PSU, SSU, FWT, URUT, WI1, WI2, STRATA,
    R703_A,
    # yang baru dimutate
    perdesaan, status_rumah, art_banyak, perempuan, usia_tua, 
    pendidikan_tinggi,pertanian,jam_kerja,keluhan,
     dapat_bpnt, dapat_blt, dapat_bp,dapat_pktd,miskin,
    fies1, fies2, fies3, fies4, fies5, fies6, fies7, fies8,
    kab_kota)

1.6 kategori rawan pangan berdasarkan fies

buang unit analisis yang menjawab kode 8 atau 9 kategorikan individu tidak rawan pangan jika menjawab semua pertanyaan poin 5 (tidak) menjadi 0

data1 %>%
  # 1. Membuang responden dengan kode 8 atau 9 di variabel FIES
  filter(across(starts_with("fies"), ~ !.x %in% c(8, 9))) %>% summarise(jumlah_obervasi = sum(FWT))

data2 <- data1 %>%
  filter(across(starts_with("fies"), ~ !.x %in% c(8, 9)))

data2 %>%
  # 2. Membuat klasifikasi ketahanan pangan
  mutate(
    status_pangan = ifelse(
      rowSums(across(starts_with("fies"), ~ .x == 5)) == 8,  # Semua jawab 5 (jumlah item FIES = 8)
      0L,  # Tahan pangan
      1L   # Selain itu: Rawan pangan
    ),
    status_pangan = factor(
      status_pangan,
      levels = c(0, 1),
      labels = c("Tahan Pangan", "Rawan Pangan")
    )
  ) %>%  select(!starts_with("fies")) %>% summarise(jumlah_obervasi = sum(FWT))

data3 <- data2 %>%
  # 2. Membuat klasifikasi ketahanan pangan
  mutate(
    status_pangan = ifelse(
      rowSums(across(starts_with("fies"), ~ .x == 5)) == 8,  # Semua jawab 5 (jumlah item FIES = 8)
      0L,  # Tahan pangan
      1L   # Selain itu: Rawan pangan
    ),
    status_pangan = factor(
      status_pangan,
      levels = c(0, 1),
      labels = c("Tahan Pangan", "Rawan Pangan")
    )
  )

nrow(data2)

1.7 data untuk level kabupaten dan kota

masukkan data lalu cek matriks korelasi untuk melihat multikolinearitas

# Membuat data frame
data_level2 <- data.frame(
  id = c(1:12, 71:77),
   panjang = c(9.22, 80.8, 98.82, 40.71, 64.05, 72.22, 114.6, 115.9, 116.71, 
          40.5, 93.68, 105.53, 22.19, 1, 22.55, 0, 0, 8.08, 1.84),
  pasar = c(4, 42, 43, 51, 41, 30, 40, 67, 41, 24, 31, 43, 10, 2, 6, 7, 5, 4, 6),
  beras = c(272.09, 97911.44, 95940.43, 33880.74, 90280.26, 78692.15, 83368.35, 
            77211.17, 75245.95, 28961.56, 27099.05, 26919.92, 28477.90, 7034.67, 
            6444.34, 2575.77, 2185.11, 13762.48, 9162.34),
  jagung = c(16.84, 127627.90, 4968.99, 5731.35, 17988.33, 36924.03, 87088.66, 
             46367.72, 96001.23, 59828.98, 4072.68, 241815.90, 49.48, 103.90, 
             1385.39, 0.00, 333.08, 2115.76, 878.13)
)

data_level2$berasjagung <- data_level2$beras+data_level2$jagung
data_level2 <- data_level2[, !(names(data_level2) %in% c("beras", "jagung"))]

cor(data_level2[,-1])

data_finale <- inner_join(data3,data_level2,by="kab_kota")
nrow(data_finale)

sum(is.na(data_finale))

1.8 design kompleks

design <- svydesign(
  id = ~PSU + SSU,          # PSU dan SSU
  strata = ~STRATA,         # Stratifikasi
  weights = ~FWT,           # Final Weight
  data = data_finale,
  nest = TRUE               # penting untuk desain bertingkat
)

2 statistik deskriptif

termasuk proporsi rawan pangan/tidak , tabel kontingensi, dan box plot untuk jam kerja

# status pangan dan keluhan
svymean(~factor(status_pangan), design)
svymean(~factor(keluhan), design)
table(data_finale$status_pangan)

# Daftar variabel prediktor kategorik
variabel_prediktor <- c(
  "perdesaan", "status_rumah", "art_banyak", "perempuan", "usia_tua",
  "pendidikan_tinggi", "pertanian", "miskin", "dapat_bpnt", "dapat_blt",
  "dapat_pktd", "dapat_bp", "keluhan"
)

# Buat list hasil proporsi
hasil_proporsi <- lapply(variabel_prediktor, function(var) {
  # Buat formula dinamis
  formula_svy <- as.formula(paste("~status_pangan +", var))
  
  # Buat tabel kontingensi dengan svytable
  tab <- svytable(formula_svy, design)
  
  # Hitung proporsi kolom
  prop <- prop.table(tab, margin = 2)
  
  list(
    variabel = var,
    tabel = prop
  )
})

for (hasil in hasil_proporsi) {
  cat("\n===== Variabel:", hasil$variabel, "=====\n")
  print(hasil$tabel)
}

## untuk keluhan 
prop.table(svytable(~keluhan + status_pangan, design), margin = 2)  # proporsi kolom

# boxplot jumlah jam kerja vs status pangan dan pertanian

# Gabungkan status_pangan dan pertanian menjadi satu variabel dengan nama gabungan
design$variables$group <- interaction(design$variables$status_pangan, 
                                      design$variables$pertanian, 
                                      sep = " - ")

# Pastikan jadi faktor agar levelnya teratur (opsional: atur urutan level manual jika mau)
design$variables$group <- factor(design$variables$group)

# Buat boxplot dengan label sumbu X yang rapi
svyboxplot(jam_kerja ~ group, 
           design = subset(design, R703_A == "A"),
           main = "",
           xlab = "Status Pangan dan Status Pertanian", 
           ylab = "Jumlah Jam Kerja per Minggu",
           col = c("forestgreen", "forestgreen", "royalblue", "royalblue")[1:length(levels(design$variables$group))])

3 regresi logistik pengaruh

3.1 pengujian random effect

alpha <- 0.1
# model logistik biner : intercept only
fit_bio <- glm(formula=status_pangan~1,family=binomial("logit"),data=data_finale)
summary(fit_bio)

# model logistik biner multilevel : intercept only
fit_mio <- glmer(formula=status_pangan~(1|kab_kota),family=binomial("logit"),data=data_finale)
summary(fit_bio)

#uji signifikansi random effect
-2*(logLik(fit_bio)-logLik(fit_mio))

#wilayah kritis
pchisq(alpha,df=1,lower.tail = F)

lrtest(fit_bio,fit_mio)

3.2 icc

random effect di multilevel intercept variancenya dibagi random effect di multilevel intercept variancenya + varians level 1 dengan hasil 3,29

(0.2065)/(0.2065+3.29)

3.3 simultan dan parsial

# model logistik biner : conditional model
fit_b <- glm(formula =status_pangan~perdesaan+perempuan+usia_tua+art_banyak+
            pendidikan_tinggi+status_rumah+pertanian+miskin+
            dapat_pktd+dapat_bp,family=binomial("logit"),
            data=data_finale)

summary(fit_b)

#uji signifikansi secara parsial : wald test
fit_mb1 <- glmer(formula =status_pangan~perdesaan+perempuan+usia_tua+art_banyak+
            pendidikan_tinggi+status_rumah+pertanian+miskin+
            dapat_pktd+dapat_bp+berasjagung+
            panjang+pasar+(1|kab_kota),family=binomial,
            control=glmerControl(optimizer = "bobyqa"),nAGQ=0,data=data_finale)

summary(fit_mb1)

lrtest(fit_mio,fit_mb1)

3.4 odds ratio

# Odds Ratio
exp_in <- as.vector(exp(summary(fit_mb1)$coefficients[1:14, 1]))
exp_in <- round(exp_in, 4)

ID <- as.vector(c("intercept", "perdesaan", "perempuan", "usia_tua", 
                  "art_banyak", "pendidikan_tinggi","status_rumah", 
                  "pertanian","miskin","PKTD","pangan/beras",
                  "berasjagung","panjang","pasar"))

odds_in <- cbind(ID, exp_in)
odds_in

3.5 kebaikan model AUC-ROC

untuk melihat seberapa baik kemampuan prediksi antara 2 model

# Model GLM biasa
prob_fit_b <- predict(fit_b, type = "response")

# Model GLMM (menggunakan glmer)
prob_fit_mb1 <- predict(fit_mb1, type = "response")

# Buat ROC dan AUC untuk model 1
roc_b <- roc(data_finale$status_pangan, prob_fit_b)
auc_b <- auc(roc_b)

# Buat ROC dan AUC untuk model 2
roc_mb1 <- roc(data_finale$status_pangan, prob_fit_mb1)
auc_mb1 <- auc(roc_mb1)

# Tampilkan hasil
print(auc_b)
print(auc_mb1)

# Ekstrak koordinat (FPR dan TPR)
roc_df_b <- data.frame(
  FPR = 1 - roc_b$specificities,
  TPR = roc_b$sensitivities,
  Model = "GLM"
)

roc_df_mb1 <- data.frame(
  FPR = 1 - roc_mb1$specificities,
  TPR = roc_mb1$sensitivities,
  Model = "GLMM"
)

# Gabungkan
roc_df_all <- rbind(roc_df_b, roc_df_mb1)

# Plot ROC + garis diagonal
ggplot(roc_df_all, aes(x = FPR, y = TPR, color = Model)) +
  geom_line(size = 1.2) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +  # garis AUC = 0.5
  scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +
  labs(title = "",
       x = "False Positive Rate (1 - Specificity)",
       y = "True Positive Rate (Sensitivity)") +
  theme_minimal() +
 theme(
  panel.grid = element_blank(),              # hapus garis kisi
  panel.background = element_blank(),        # hapus latar panel
  plot.background = element_blank(),         # hapus latar plot
  axis.line = element_line(color = "gray")  # tampilkan garis sumbu x & y
  ) +
  guides(color = guide_legend(title = NULL))

4 untuk dampak

4.1 anova untuk jam kerja

## cek centroid weight
# Rata-rata per kelompok status_pangan dan sektor, menggunakan bobot survei
svyby(~jam_kerja, ~status_pangan + pertanian, design, svymean)

# untuk robust
data_finale %>%
  filter(R703_A == "A") %>%   
  group_by(pertanian, status_pangan) %>%
  summarise(
    rata_trimmed = mean(jam_kerja, trim = 0.2),
    rata_biasa = mean(jam_kerja),
    .groups = "drop"
  )

# Robust ANOVA (trimmed mean 20% default)
t1way(jam_kerja ~ status_pangan, data = data_finale[data_finale$pertanian=="Pertanian",], tr = 0.2)
t1way(jam_kerja ~ status_pangan, data = data_finale[data_finale$pertanian=="NonPertanian",], tr = 0.2)

4.2 logistik biner

4.2.1 parsial

##  model intercept 
fit_in <- svyglm(keluhan ~ 1, design = design, family = quasibinomial())

## modul full
fit_ls <- svyglm(keluhan ~ status_pangan, design = design, family = quasibinomial())

summary(fit_ls)

4.2.2 odds ratio

# Odds Ratio
exp_1 <- as.vector(exp(summary(fit_ls)$coefficients[1:2, 1]))
exp_1 <- round(exp_1, 4)

ID <- as.vector(c("intercept", "rawan"))

odds_1 <- cbind(ID, exp_1)
odds_1