Pengolahan Setelah SE

1 Load Library

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(haven)
library(purrr)
library(stringr)
library(tibble)
library(knitr)
library(tidyr)
library(forcats)
library(ggplot2)
library(openxlsx)

## Warning: package 'openxlsx' was built under R version 4.3.3

library(logistf)

## Warning: package 'logistf' was built under R version 4.3.3

## Registered S3 method overwritten by 'formula.tools':
##   method               from    
##   as.character.formula openxlsx

2 Import Dataset

library(haven)
library(dplyr)

path_data_join_all <- "C:/Users/nisri/Documents/KULIAH/SEMESTER 7/ZKRIPZWIT/POSTPARTUM DEPRESSION/OLAH DATA/After Sempro/(V1) data_join_all_hasil bismillah V1.sav" 

data_all <- read_sav(path_data_join_all)

# Cek apakah sudah jadi data frame
class(data_all)

## [1] "tbl_df"     "tbl"        "data.frame"

dim(data_all)

## [1] 64034   121

View(data_all)

3 Filter Unit Analisis

3.1 Provinsi Jawa Barat

# ============================================================
# FILTER UNIT ANALISIS
# ============================================================

## Provinsi Jawa Barat dan Jawa Tengah

# Cek jumlah data awal
dim(data_all)

## [1] 64034   121

# Cek isi kode provinsi
data_all %>%
  count(B1R1, sort = TRUE)

# Filter data Provinsi Jawa Barat dan Jawa Tengah
# 32 = Jawa Barat
# 33 = Jawa Tengah

data_jabar <- data_all %>%
  filter(as.numeric(B1R1) %in% c(32))

# Cek hasil filter provinsi
dim(data_jabar)

## [1] 4281  121

data_jabar %>%
  count(B1R1, sort = TRUE)

3.2 Wilayah Perkotaan pada Provinsi Jawa Barat

## Wilayah Perkotaan pada Provinsi Jawa Barat dan Jawa Tengah

# Cek isi kode wilayah
data_jabar %>%
  count(B1R5, sort = TRUE)

# Filter data Jabar-Jateng khusus perkotaan
# 1 = Perkotaan
# 2 = Perdesaan

data_jabar_kota <- data_jabar %>%
  filter(as.numeric(B1R5) == 1)

# Cek hasil filter Jabar-Jateng perkotaan
dim(data_jabar_kota)

## [1] 3619  121

data_jabar_kota %>%
  count(B1R1, B1R5, sort = TRUE)

data_all

# data awal hasil import

data_jabar

# data_all yang barisnya hanya Provinsi Jawa Barat 

data_jabar_kota

# data_jabar yang barisnya hanya wilayah perkotaan

3.3 Ibu dengan Balita <= 12 bulan

# ============================================================
# FILTER UNIT ANALISIS: ANAK/BALITA USIA <= 12 BULAN
# ============================================================

# Cek dulu nama variabel yang berkaitan dengan tanggal/umur
names(data_jabar_kota)[
  grepl("tgl|tanggal|bln|bulan|thn|tahun|umur|lahir", 
        names(data_jabar_kota), 
        ignore.case = TRUE)
]

##  [1] "B4K7BLN"                 "B4K7THN"                
##  [3] "B4K7BLN_balita"          "tgl_lahir_balita_8digit"
##  [5] "tgl_wawancara_8digit"    "tgl_lahir_balita"       
##  [7] "tgl_wawancara"           "umur_balita_hari"       
##  [9] "umur_balita_bulan"       "umur_0_6_bulan"         
## [11] "umur_0_12_bulan"         "umur_0_24_bulan"

B4K6 dan B2R2 digunakan untuk mengetahui umur balita (bulan dan hari)

# ============================================================
# FILTER UNIT ANALISIS: IBU DENGAN BALITA <= 12 BULAN
# ============================================================

library(stringr)
library(lubridate)

## Warning: package 'lubridate' was built under R version 4.3.2

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

# Fungsi singkat untuk ubah tanggal SKI ke Date
ubah_tanggal_ski <- function(x) {
  x <- haven::zap_labels(x)
  x <- as.character(x)
  x <- str_replace(x, "\\.0$", "")
  x <- str_pad(x, width = 8, side = "left", pad = "0")
  as.Date(x, format = "%d%m%Y")
}

# Membuat variabel bantu tanggal dan umur dari variabel mentah
data_jabar_kota_umur <- data_jabar_kota %>%
  mutate(
    tgl_lahir_balita_date = ubah_tanggal_ski(B4K6),
    tgl_wawancara_date    = ubah_tanggal_ski(B2R2),
    
    umur_balita_hari = as.integer(tgl_wawancara_date - tgl_lahir_balita_date),
    
    umur_balita_bulan = 
      (year(tgl_wawancara_date) - year(tgl_lahir_balita_date)) * 12 +
      (month(tgl_wawancara_date) - month(tgl_lahir_balita_date)) -
      ifelse(day(tgl_wawancara_date) < day(tgl_lahir_balita_date), 1, 0)
  )

# Cek hasil hitung umur
data_jabar_kota_umur %>%
  dplyr::select(
    B4K6,
    B2R2,
    tgl_lahir_balita_date,
    tgl_wawancara_date,
    umur_balita_hari,
    umur_balita_bulan
  ) %>%
  head(30)

view(data_jabar_kota_umur)

# Cek ringkasan umur
data_jabar_kota_umur %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    umur_missing = sum(is.na(umur_balita_bulan)),
    umur_negatif = sum(umur_balita_bulan < 0, na.rm = TRUE),
    umur_min = min(umur_balita_bulan, na.rm = TRUE),
    umur_max = max(umur_balita_bulan, na.rm = TRUE)
  )

filter <=12 bulan

# Filter ibu dengan balita usia <= 12 bulan
data_jabar_kota_12 <- data_jabar_kota_umur %>%
  dplyr::filter(
    !is.na(umur_balita_bulan),
    umur_balita_bulan >= 0,
    umur_balita_bulan <= 12
  )

# Cek hasil filter
dim(data_jabar_kota)

## [1] 3619  121

dim(data_jabar_kota_12)

## [1] 734 123

# Distribusi umur balita setelah filter
data_jabar_kota_12 %>%
  dplyr::count(umur_balita_bulan, sort = FALSE)

# Cek provinsi dan wilayah
data_jabar_kota_12 %>%
  dplyr::count(B1R1, B1R5, sort = TRUE)

View(data_jabar_kota_12)

data_all → data_jabar → data_jabar_kota → data_jabar_kota_umur → data_jabar_kota_12

4 Variabel Y: Status Depresi

# ============================================================
# CEK MISSING VARIABEL DEPRESI
# ============================================================

data_jabar_kota_12 %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_JMLC01_03 = sum(is.na(JMLC01_03_depresi)),
    missing_JMLC04_10 = sum(is.na(JMLC04_10_depresi)),
    missing_salah_satu = sum(
      is.na(JMLC01_03_depresi) | is.na(JMLC04_10_depresi)
    ),
    lengkap_keduanya = sum(
      !is.na(JMLC01_03_depresi) & !is.na(JMLC04_10_depresi)
    )
  )

# status_depresi dibuat berdasarkan kriteria minimal
# 2 gejala pada JMLC01_03_depresi dan minimal 2 gejala pada JMLC04_10_depresi.
# ============================================================
# MEMBUAT DATA DEPRESI DAN VARIABEL Y
# ============================================================

data_depresi <- data_jabar_kota_12 %>%
  dplyr::filter(
    !is.na(JMLC01_03_depresi),
    !is.na(JMLC04_10_depresi)
  ) %>%
  dplyr::mutate(
    status_depresi = dplyr::case_when(
      JMLC01_03_depresi >= 2 & JMLC04_10_depresi >= 2 ~ 1,
      TRUE ~ 0
    )
  )

# ============================================================
# CEK HASIL DATA DEPRESI
# ============================================================

dim(data_jabar_kota_12)

## [1] 734 123

dim(data_depresi)

## [1] 733 124

# Pastikan sudah tidak ada missing pada variabel depresi
data_depresi %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_JMLC01_03 = sum(is.na(JMLC01_03_depresi)),
    missing_JMLC04_10 = sum(is.na(JMLC04_10_depresi)),
    missing_status_depresi = sum(is.na(status_depresi))
  )

# Cek distribusi status depresi
data_depresi %>%
  dplyr::count(status_depresi) %>%
  dplyr::mutate(
    persen = round(n / sum(n) * 100, 2)
  )

# ============================================================
# MEMBUAT DATA ANALISIS
# ============================================================

data_analisis <- data_depresi

5 Variabel X

5.1 Variabel X1: H02B (Paritas)

# ============================================================
# MEMBUAT X1: PARITAS
# ============================================================

data_analisis <- data_depresi %>%
  dplyr::mutate(
    H02B_num = as.numeric(haven::zap_labels(H02B)),
    
    X1 = dplyr::case_when(
      is.na(H02B_num) ~ NA_real_,
      H02B_num >= 2 ~ 0,
      H02B_num < 2 ~ 1
    )
  )

# Cek distribusi X1
data_analisis %>%
  dplyr::count(X1) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

5.2 Variabel X2: H32 (Metode Persalinan)

# ============================================================
# X2: METODE PERSALINAN
# ============================================================
# H32:
# 1 = Persalinan normal
# 2, 3, 4, 5 = Selain normal
#
# X2:
# 0 = Normal
# 1 = Selain normal
# ============================================================

# Cek distribusi variabel asli H32
data_analisis %>%
  dplyr::mutate(
    H32_num = as.numeric(haven::zap_labels(H32))
  ) %>%
  dplyr::count(H32, H32_num, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 4 × 3
##   H32               H32_num     n
##   <dbl+lbl>           <dbl> <int>
## 1 1 [Normal]              1   525
## 2 2 [Operasi sesar]       2   202
## 3 3 [Vacuum]              3     5
## 4 5 [Lainnya]             5     1

# Cek missing dan kode khusus H32
data_analisis %>%
  dplyr::mutate(
    H32_num = as.numeric(haven::zap_labels(H32))
  ) %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_H32 = sum(is.na(H32_num)),
    kode_1_normal = sum(H32_num == 1, na.rm = TRUE),
    kode_2_5_selain_normal = sum(H32_num %in% c(2, 3, 4, 5), na.rm = TRUE),
    kode_lain = sum(!is.na(H32_num) & !(H32_num %in% c(1, 2, 3, 4, 5)))
  )

# Membuat X2 pada data_analisis
data_analisis <- data_analisis %>%
  dplyr::mutate(
    H32_num = as.numeric(haven::zap_labels(H32)),
    
    X2 = dplyr::case_when(
      is.na(H32_num) ~ NA_real_,
      H32_num == 1 ~ 0,
      H32_num %in% c(2, 3, 4, 5) ~ 1,
      TRUE ~ NA_real_
    )
  )

# Label untuk SPSS
data_analisis$X2 <- haven::labelled(
  data_analisis$X2,
  labels = c(
    "Normal" = 0,
    "Selain normal" = 1
  ),
  label = "Metode persalinan"
)

# Cek recoding H32 menjadi X2
data_analisis %>%
  dplyr::count(H32, H32_num, X2, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 4 × 4
##   H32               H32_num X2                    n
##   <dbl+lbl>           <dbl> <dbl+lbl>         <int>
## 1 1 [Normal]              1 0 [Normal]          525
## 2 2 [Operasi sesar]       2 1 [Selain normal]   202
## 3 3 [Vacuum]              3 1 [Selain normal]     5
## 4 5 [Lainnya]             5 1 [Selain normal]     1

# Cek distribusi X2
data_analisis %>%
  dplyr::count(X2) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

5.3 Variabel X3: H46 (Komplikasi Masa Nifas)

# ============================================================
# X3: KOMPLIKASI NIFAS
# ============================================================
# X3:
# 0 = Tidak ada komplikasi nifas
# 1 = Ada komplikasi nifas
# ============================================================


# CEK ISI VARIABEL H46
# Cek distribusi H46 apa adanya
data_analisis %>%
  dplyr::count(H46, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 31 × 2
##    H46               n
##    <chr>         <int>
##  1 "          Z"   628
##  2 "      G"        19
##  3 "       H"       13
##  4 "        I"      12
##  5 "A"               9
##  6 "   D"            8
##  7 "         J"      6
##  8 "     F"          4
##  9 "  C"             4
## 10 "      GH"        3
## 11 "     FG"         3
## 12 " B"              3
## 13 "  CD"            2
## 14 "A  D"            2
## 15 "      G I"       1
## 16 "     FGH"        1
## 17 "   D  G"         1
## 18 "   D  GH"        1
## 19 "   D F H"        1
## 20 "   D FG"         1
## 21 "   D FG  J"      1
## 22 "  C     I"       1
## 23 "  C    HI"       1
## 24 "  C   GH"        1
## 25 "  CD   H"        1
## 26 " BCD  G"         1
## 27 "A  D   H"        1
## 28 "A  D F HI"       1
## 29 "A CDEFGH"        1
## 30 "AB     H"        1
## 31 "AB D    I"       1

# Cek missing H46
data_analisis %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_H46 = sum(is.na(H46))
  )

Karena data nya banyak spasi, bikin versi bersih dahulu agar mudah di recode

# ============================================================
# X3: KOMPLIKASI NIFAS
# ============================================================
# H46:
# Z = Tidak ada komplikasi nifas
# A-J = Ada komplikasi nifas
#
# X3:
# 0 = Tidak ada komplikasi nifas
# 1 = Ada komplikasi nifas
# ============================================================

data_analisis <- data_analisis %>%
  dplyr::mutate(
    H46_clean = stringr::str_replace_all(H46, "\\s+", ""),
    
    X3 = dplyr::case_when(
      is.na(H46_clean) | H46_clean == "" ~ NA_real_,
      H46_clean == "Z" ~ 0,
      stringr::str_detect(H46_clean, "[A-J]") ~ 1,
      TRUE ~ NA_real_
    )
  )

# Cek recoding H46 menjadi X3
data_analisis %>%
  dplyr::count(H46, H46_clean, X3, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 31 × 4
##    H46           H46_clean    X3     n
##    <chr>         <chr>     <dbl> <int>
##  1 "          Z" Z             0   628
##  2 "      G"     G             1    19
##  3 "       H"    H             1    13
##  4 "        I"   I             1    12
##  5 "A"           A             1     9
##  6 "   D"        D             1     8
##  7 "         J"  J             1     6
##  8 "     F"      F             1     4
##  9 "  C"         C             1     4
## 10 "      GH"    GH            1     3
## 11 "     FG"     FG            1     3
## 12 " B"          B             1     3
## 13 "  CD"        CD            1     2
## 14 "A  D"        AD            1     2
## 15 "      G I"   GI            1     1
## 16 "     FGH"    FGH           1     1
## 17 "   D  G"     DG            1     1
## 18 "   D  GH"    DGH           1     1
## 19 "   D F H"    DFH           1     1
## 20 "   D FG"     DFG           1     1
## 21 "   D FG  J"  DFGJ          1     1
## 22 "  C     I"   CI            1     1
## 23 "  C    HI"   CHI           1     1
## 24 "  C   GH"    CGH           1     1
## 25 "  CD   H"    CDH           1     1
## 26 " BCD  G"     BCDG          1     1
## 27 "A  D   H"    ADH           1     1
## 28 "A  D F HI"   ADFHI         1     1
## 29 "A CDEFGH"    ACDEFGH       1     1
## 30 "AB     H"    ABH           1     1
## 31 "AB D    I"   ABDI          1     1

# Cek distribusi X3
data_analisis %>%
  dplyr::count(X3) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Label untuk SPSS
data_analisis$X3 <- haven::labelled(
  data_analisis$X3,
  labels = c(
    "Tidak ada komplikasi nifas" = 0,
    "Ada komplikasi nifas" = 1
  ),
  label = "Komplikasi nifas"
)

5.4 Variabel X4: H04E (Status Kehamilan)

# ============================================================
# X4: STATUS KEHAMILAN DIINGINKAN
# ============================================================
# H04E:
# 1 = Diinginkan saat itu
# 2 = Diinginkan kemudian
# 3 = Tidak diinginkan
#
# X4:
# 0 = Diinginkan saat itu
# 1 = Diinginkan kemudian / tidak diinginkan
# ============================================================

# Cek distribusi variabel asli H04E
data_analisis %>%
  dplyr::mutate(
    H04E_num = as.numeric(haven::zap_labels(H04E)),
    H04E_label = as.character(haven::as_factor(H04E))
  ) %>%
  dplyr::count(H04E, H04E_num, H04E_label, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 3 × 4
##   H04E                    H04E_num H04E_label              n
##   <dbl+lbl>                  <dbl> <chr>               <int>
## 1 1 [Diinginkan saat itu]        1 Diinginkan saat itu   620
## 2 2 [Diinginkan kemudian]        2 Diinginkan kemudian    86
## 3 3 [Tidak diinginkan]           3 Tidak diinginkan       27

# Cek missing dan kode khusus H04E
data_analisis %>%
  dplyr::mutate(
    H04E_num = as.numeric(haven::zap_labels(H04E))
  ) %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_H04E = sum(is.na(H04E_num)),
    kode_1_diinginkan_saat_itu = sum(H04E_num == 1, na.rm = TRUE),
    kode_2_3_tidak_saat_itu = sum(H04E_num %in% c(2, 3), na.rm = TRUE),
    kode_lain = sum(!is.na(H04E_num) & !(H04E_num %in% c(1, 2, 3)))
  )

# Membuat X4 pada data_analisis
data_analisis <- data_analisis %>%
  dplyr::mutate(
    H04E_num = as.numeric(haven::zap_labels(H04E)),
    
    X4 = dplyr::case_when(
      H04E_num == 1 ~ 0,
      H04E_num %in% c(2, 3) ~ 1,
      TRUE ~ NA_real_
    )
  )

# Cek recoding H04E menjadi X4
data_analisis %>%
  dplyr::count(H04E, H04E_num, X4, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 3 × 4
##   H04E                    H04E_num    X4     n
##   <dbl+lbl>                  <dbl> <dbl> <int>
## 1 1 [Diinginkan saat itu]        1     0   620
## 2 2 [Diinginkan kemudian]        2     1    86
## 3 3 [Tidak diinginkan]           3     1    27

# Cek distribusi X4
data_analisis %>%
  dplyr::count(X4) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Label untuk SPSS
data_analisis$X4 <- haven::labelled(
  data_analisis$X4,
  labels = c(
    "Diinginkan saat itu" = 0,
    "Diinginkan kemudian/tidak diinginkan" = 1
  ),
  label = "Status kehamilan diinginkan"
)

5.5 Variabel X5: B4K7THN (Umur Ibu)

# ============================================================
# X5: USIA IBU
# ============================================================
# B4K7THN = umur ibu
#
# X5:
# 0 = Umur ibu 20-34 tahun
# 1 = Umur ibu <20 tahun atau >=35 tahun
# ============================================================

# Cek distribusi variabel asli B4K7THN
data_analisis %>%
  dplyr::mutate(
    B4K7THN_num = as.numeric(haven::zap_labels(B4K7THN))
  ) %>%
  dplyr::count(B4K7THN, B4K7THN_num, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 31 × 3
##    B4K7THN   B4K7THN_num     n
##    <dbl+lbl>       <dbl> <int>
##  1 30                 30    55
##  2 32                 32    54
##  3 28                 28    52
##  4 31                 31    45
##  5 26                 26    43
##  6 33                 33    39
##  7 23                 23    38
##  8 27                 27    37
##  9 35                 35    36
## 10 38                 38    34
## 11 29                 29    33
## 12 25                 25    31
## 13 34                 34    31
## 14 39                 39    30
## 15 22                 22    24
## 16 37                 37    24
## 17 36                 36    23
## 18 24                 24    20
## 19 21                 21    19
## 20 40                 40    15
## 21 41                 41     9
## 22 43                 43     9
## 23 44                 44     8
## 24 19                 19     5
## 25 42                 42     5
## 26 18                 18     4
## 27 20                 20     4
## 28 16                 16     2
## 29 45                 45     2
## 30 46                 46     1
## 31 47                 47     1

# Cek missing dan kelompok umur ibu
data_analisis %>%
  dplyr::mutate(
    B4K7THN_num = as.numeric(haven::zap_labels(B4K7THN))
  ) %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_B4K7THN = sum(is.na(B4K7THN_num)),
    umur_min = min(B4K7THN_num, na.rm = TRUE),
    umur_max = max(B4K7THN_num, na.rm = TRUE),
    umur_kurang_20 = sum(B4K7THN_num < 20, na.rm = TRUE),
    umur_20_34 = sum(B4K7THN_num >= 20 & B4K7THN_num <= 34, na.rm = TRUE),
    umur_35_keatas = sum(B4K7THN_num >= 35, na.rm = TRUE)
  )

# Membuat X5 pada data_analisis
data_analisis <- data_analisis %>%
  dplyr::mutate(
    B4K7THN_num = as.numeric(haven::zap_labels(B4K7THN)),
    
    X5 = dplyr::case_when(
      B4K7THN_num >= 20 & B4K7THN_num <= 34 ~ 0,
      B4K7THN_num < 20 | B4K7THN_num >= 35 ~ 1,
      TRUE ~ NA_real_
    )
  )

# Cek recoding B4K7THN menjadi X5
data_analisis %>%
  dplyr::count(B4K7THN, B4K7THN_num, X5, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 31 × 4
##    B4K7THN   B4K7THN_num    X5     n
##    <dbl+lbl>       <dbl> <dbl> <int>
##  1 30                 30     0    55
##  2 32                 32     0    54
##  3 28                 28     0    52
##  4 31                 31     0    45
##  5 26                 26     0    43
##  6 33                 33     0    39
##  7 23                 23     0    38
##  8 27                 27     0    37
##  9 35                 35     1    36
## 10 38                 38     1    34
## 11 29                 29     0    33
## 12 25                 25     0    31
## 13 34                 34     0    31
## 14 39                 39     1    30
## 15 22                 22     0    24
## 16 37                 37     1    24
## 17 36                 36     1    23
## 18 24                 24     0    20
## 19 21                 21     0    19
## 20 40                 40     1    15
## 21 41                 41     1     9
## 22 43                 43     1     9
## 23 44                 44     1     8
## 24 19                 19     1     5
## 25 42                 42     1     5
## 26 18                 18     1     4
## 27 20                 20     0     4
## 28 16                 16     1     2
## 29 45                 45     1     2
## 30 46                 46     1     1
## 31 47                 47     1     1

# Cek distribusi X5
data_analisis %>%
  dplyr::count(X5) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Label untuk SPSS
data_analisis$X5 <- haven::labelled(
  data_analisis$X5,
  labels = c(
    "20-34 tahun" = 0,
    "<20 tahun atau >=35 tahun" = 1
  ),
  label = "Umur ibu"
)

5.5.1 VAR X5_Num

# Membuat X5_num pada data_analisis
data_analisis <- data_analisis %>%
  dplyr::mutate(
    B4K7THN_num = as.numeric(haven::zap_labels(B4K7THN)),
    
    X5_num = B4K7THN_num
  )

5.6 Variabel X6: B4K8 (Pendidikan Ibu)

# ============================================================
# X6: PENDIDIKAN IBU
# ============================================================
# B4K8 = pendidikan ibu
#
# X6:
# 0 = SLTA atau lebih tinggi
# 1 = SMP atau lebih rendah
# ============================================================

# Cek distribusi variabel asli B4K8
data_analisis %>%
  dplyr::mutate(
    B4K8_num = as.numeric(haven::zap_labels(B4K8)),
    B4K8_label = as.character(haven::as_factor(B4K8))
  ) %>%
  dplyr::count(B4K8, B4K8_num, B4K8_label, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 7 × 4
##   B4K8                            B4K8_num B4K8_label                      n
##   <dbl+lbl>                          <dbl> <chr>                       <int>
## 1 5 [Tamat SLTA/MA]                      5 Tamat SLTA/MA                 322
## 2 4 [Tamat SLTP/MTS]                     4 Tamat SLTP/MTS                177
## 3 3 [Tamat SD/MI]                        3 Tamat SD/MI                    89
## 4 7 [Tamat PT]                           7 Tamat PT                       75
## 5 6 [Tamat D1/D2/D3]                     6 Tamat D1/D2/D3                 43
## 6 2 [Tidak tamat SD/MI]                  2 Tidak tamat SD/MI              25
## 7 1 [Tidak/ belum pernah sekolah]        1 Tidak/ belum pernah sekolah     2

# Cek missing dan kode pendidikan ibu
data_analisis %>%
  dplyr::mutate(
    B4K8_num = as.numeric(haven::zap_labels(B4K8))
  ) %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_B4K8 = sum(is.na(B4K8_num)),
    SMP_atau_lebih_rendah = sum(B4K8_num %in% c(1, 2, 3, 4), na.rm = TRUE),
    SLTA_atau_lebih_tinggi = sum(B4K8_num %in% c(5, 6, 7), na.rm = TRUE),
    kode_lain = sum(!is.na(B4K8_num) & !(B4K8_num %in% c(1, 2, 3, 4, 5, 6, 7)))
  )

# Membuat X6 pada data_analisis
data_analisis <- data_analisis %>%
  dplyr::mutate(
    B4K8_num = as.numeric(haven::zap_labels(B4K8)),
    
    X6 = dplyr::case_when(
      B4K8_num %in% c(5, 6, 7) ~ 0,
      B4K8_num %in% c(1, 2, 3, 4) ~ 1,
      TRUE ~ NA_real_
    )
  )

# Cek recoding B4K8 menjadi X6
data_analisis %>%
  dplyr::count(B4K8, B4K8_num, X6, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 7 × 4
##   B4K8                            B4K8_num    X6     n
##   <dbl+lbl>                          <dbl> <dbl> <int>
## 1 5 [Tamat SLTA/MA]                      5     0   322
## 2 4 [Tamat SLTP/MTS]                     4     1   177
## 3 3 [Tamat SD/MI]                        3     1    89
## 4 7 [Tamat PT]                           7     0    75
## 5 6 [Tamat D1/D2/D3]                     6     0    43
## 6 2 [Tidak tamat SD/MI]                  2     1    25
## 7 1 [Tidak/ belum pernah sekolah]        1     1     2

# Cek distribusi X6
data_analisis %>%
  dplyr::count(X6) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Label untuk SPSS
data_analisis$X6 <- haven::labelled(
  data_analisis$X6,
  labels = c(
    "SLTA atau lebih tinggi" = 0,
    "SMP atau lebih rendah" = 1
  ),
  label = "Pendidikan ibu"
)

5.7 Variabel X7: Kuintil (Status Ekonomi)

# ============================================================
# X7: STATUS EKONOMI
# ============================================================
# Kuintil:
# 1, 2 = Miskin
# 3, 4, 5 = Menengah-kaya
#
# X9:
# 0 = Menengah-kaya (Q3-Q5)
# 1 = Miskin (Q1-Q2)
# ============================================================

# Cek distribusi variabel asli Kuintil
data_analisis %>%
  dplyr::mutate(
    Kuintil_num = as.numeric(haven::zap_labels(Kuintil))
  ) %>%
  dplyr::count(Kuintil, Kuintil_num, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 5 × 3
##   Kuintil            Kuintil_num     n
##   <dbl+lbl>                <dbl> <int>
## 1 5 [Teratas]                  5   269
## 2 4 [Menengah atas]            4   219
## 3 3 [Menengah]                 3   153
## 4 2 [Menengah bawah]           2    71
## 5 1 [Terbawah]                 1    21

# Cek missing dan kategori Kuintil
data_analisis %>%
  dplyr::mutate(
    Kuintil_num = as.numeric(haven::zap_labels(Kuintil))
  ) %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_Kuintil = sum(is.na(Kuintil_num)),
    miskin_Q1_Q2 = sum(Kuintil_num %in% c(1, 2), na.rm = TRUE),
    menengah_kaya_Q3_Q5 = sum(Kuintil_num %in% c(3, 4, 5), na.rm = TRUE),
    kode_lain = sum(!is.na(Kuintil_num) & !(Kuintil_num %in% c(1, 2, 3, 4, 5)))
  )

# Membuat X7 pada data_analisis
data_analisis <- data_analisis %>%
  dplyr::mutate(
    Kuintil_num = as.numeric(haven::zap_labels(Kuintil)),
    
    X7 = dplyr::case_when(
      Kuintil_num %in% c(3, 4, 5) ~ 0,
      Kuintil_num %in% c(1, 2) ~ 1,
      TRUE ~ NA_real_
    )
  )

# Cek recoding Kuintil menjadi X7
data_analisis %>%
  dplyr::count(Kuintil, Kuintil_num, X7, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 5 × 4
##   Kuintil            Kuintil_num    X7     n
##   <dbl+lbl>                <dbl> <dbl> <int>
## 1 5 [Teratas]                  5     0   269
## 2 4 [Menengah atas]            4     0   219
## 3 3 [Menengah]                 3     0   153
## 4 2 [Menengah bawah]           2     1    71
## 5 1 [Terbawah]                 1     1    21

# Cek distribusi X7
data_analisis %>%
  dplyr::count(X7) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Label untuk SPSS
data_analisis$X7 <- haven::labelled(
  data_analisis$X7,
  labels = c(
    "Menengah-kaya (Q3-Q5)" = 0,
    "Miskin (Q1-Q2)" = 1
  ),
  label = "Status ekonomi"
)

5.7.1 VAR X7_1

# Distribusi Kuintil
data_analisis %>%
  dplyr::count(Kuintil_num) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Kuintil terhadap status depresi
data_analisis %>%
  dplyr::count(Kuintil_num, status_depresi) %>%
  dplyr::group_by(Kuintil_num) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Persentase depresi per Kuintil
data_analisis %>%
  dplyr::group_by(Kuintil_num) %>%
  dplyr::summarise(
    total = dplyr::n(),
    depresi = sum(status_depresi == 1, na.rm = TRUE),
    tidak_depresi = sum(status_depresi == 0, na.rm = TRUE),
    persen_depresi = round(100 * depresi / total, 2),
    .groups = "drop"
  )

# ============================================================
# X7_1: STATUS EKONOMI ALTERNATIF
# ============================================================
# X7_1:
# 0 = Q2-Q5
# 1 = Q1
# ============================================================

data_analisis <- data_analisis %>%
  dplyr::mutate(
    X7_1 = dplyr::case_when(
      Kuintil_num %in% c(2,3,4, 5) ~ 0,
      Kuintil_num %in% c(1) ~ 1,
      TRUE ~ NA_real_
    )
  )

data_analisis %>%
  dplyr::count(X7_1) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

5.8 Variabel X8: B3R1 (Jumlah Anggota Ruta)

# ============================================================
# X8: JUMLAH ANGGOTA RUMAH TANGGA
# ============================================================
# B3R1 = jumlah anggota rumah tangga
#
# X8:
# 0 = 1-4 ART
# 1 = >=5 ART
# ============================================================

# Cek distribusi variabel asli B3R1
data_analisis %>%
  dplyr::mutate(
    B3R1_num = as.numeric(haven::zap_labels(B3R1))
  ) %>%
  dplyr::count(B3R1, B3R1_num, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 11 × 3
##     B3R1 B3R1_num     n
##    <dbl>    <dbl> <int>
##  1     4        4   281
##  2     5        5   203
##  3     3        3   159
##  4     6        6    62
##  5     7        7    15
##  6     8        8     5
##  7     2        2     4
##  8     9        9     1
##  9    10       10     1
## 10    11       11     1
## 11    13       13     1

# Cek missing dan kategori jumlah ART
data_analisis %>%
  dplyr::mutate(
    B3R1_num = as.numeric(haven::zap_labels(B3R1))
  ) %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_B3R1 = sum(is.na(B3R1_num)),
    ART_1_4 = sum(B3R1_num >= 1 & B3R1_num <= 4, na.rm = TRUE),
    ART_5_keatas = sum(B3R1_num >= 5, na.rm = TRUE),
    ART_kurang_1 = sum(B3R1_num < 1, na.rm = TRUE)
  )

# Membuat X8 pada data_analisis
data_analisis <- data_analisis %>%
  dplyr::mutate(
    B3R1_num = as.numeric(haven::zap_labels(B3R1)),
    
    X8 = dplyr::case_when(
      B3R1_num >= 1 & B3R1_num <= 4 ~ 0,
      B3R1_num >= 5 ~ 1,
      TRUE ~ NA_real_
    )
  )

# Cek recoding B3R1 menjadi X8
data_analisis %>%
  dplyr::count(B3R1, B3R1_num, X8, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 11 × 4
##     B3R1 B3R1_num    X8     n
##    <dbl>    <dbl> <dbl> <int>
##  1     4        4     0   281
##  2     5        5     1   203
##  3     3        3     0   159
##  4     6        6     1    62
##  5     7        7     1    15
##  6     8        8     1     5
##  7     2        2     0     4
##  8     9        9     1     1
##  9    10       10     1     1
## 10    11       11     1     1
## 11    13       13     1     1

# Cek distribusi X8
data_analisis %>%
  dplyr::count(X8) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Label untuk SPSS
data_analisis$X8 <- haven::labelled(
  data_analisis$X8,
  labels = c(
    "1-4 ART" = 0,
    ">=5 ART" = 1
  ),
  label = "Jumlah anggota rumah tangga"
)

5.9 Variabel X9: B3R2 (Jumlah Balita dalam Ruta)

# ============================================================
# X9: JUMLAH BALITA DALAM RUMAH TANGGA
# ============================================================
# B3R2 = jumlah balita dalam rumah tangga
#
# X9
# 0 = 1 balita
# 1 = >=2 balita
# ============================================================

# Cek distribusi variabel asli B3R2
data_analisis %>%
  dplyr::mutate(
    B3R2_num = as.numeric(haven::zap_labels(B3R2))
  ) %>%
  dplyr::count(B3R2, B3R2_num, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 4 × 3
##    B3R2 B3R2_num     n
##   <dbl>    <dbl> <int>
## 1     1        1   553
## 2     2        2   170
## 3     3        3     9
## 4     4        4     1

# Cek missing dan kategori jumlah balita
data_analisis %>%
  dplyr::mutate(
    B3R2_num = as.numeric(haven::zap_labels(B3R2))
  ) %>%
  dplyr::summarise(
    total_data = dplyr::n(),
    missing_B3R2 = sum(is.na(B3R2_num)),
    balita_1 = sum(B3R2_num == 1, na.rm = TRUE),
    balita_2_keatas = sum(B3R2_num >= 2, na.rm = TRUE),
    balita_kurang_1 = sum(B3R2_num < 1, na.rm = TRUE)
  )

# Membuat X9 pada data_analisis
data_analisis <- data_analisis %>%
  dplyr::mutate(
    B3R2_num = as.numeric(haven::zap_labels(B3R2)),
    
    X9 = dplyr::case_when(
      B3R2_num == 1 ~ 0,
      B3R2_num >= 2 ~ 1,
      TRUE ~ NA_real_
    )
  )

# Cek recoding B3R2 menjadi X9
data_analisis %>%
  dplyr::count(B3R2, B3R2_num, X9, sort = TRUE) %>%
  print(n = Inf)

## # A tibble: 4 × 4
##    B3R2 B3R2_num    X9     n
##   <dbl>    <dbl> <dbl> <int>
## 1     1        1     0   553
## 2     2        2     1   170
## 3     3        3     1     9
## 4     4        4     1     1

# Cek distribusi X9
data_analisis %>%
  dplyr::count(X9) %>%
  dplyr::mutate(
    persen = round(100 * n / sum(n), 2)
  )

# Label untuk SPSS
data_analisis$X9 <- haven::labelled(
  data_analisis$X9,
  labels = c(
    "1 balita" = 0,
    ">=2 balita" = 1
  ),
  label = "Jumlah balita dalam rumah tangga"
)

6 Analisis Deskriptif

Data awal yang digunakan

View(data_analisis)

Variabel fix yang digunakan

data_analisis <- data_analisis %>%
  dplyr::mutate(
    X1_PAR  = X1,
    X2_MP = X2,
    X3_KMN  = X3,
    X4_SKH = X4,
    X5_UI = X5,
    X5_NUM = X5_num,
    X6_PI = X6,
    X7_WI = X7,
    X8_JART = X8,
    X9_JB= X9
  )

6.1 Visualisasi Variabel Y (dengan bobot)

data_bobot <- data_analisis

persentase_ppd_bobot <- data_bobot %>%
  dplyr::filter(!is.na(status_depresi)) %>%
  dplyr::group_by(status_depresi) %>%
  dplyr::summarise(
    n_sampel = dplyr::n(),
    estimasi_populasi = sum(w_final, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  dplyr::mutate(
    status = dplyr::case_when(
      status_depresi == 0 ~ "Tidak PPD",
      status_depresi == 1 ~ "PPD"
    ),

    persen_bobot =
      100 * estimasi_populasi /
      sum(estimasi_populasi)
  ) %>%
  dplyr::mutate(
    estimasi_populasi = round(estimasi_populasi, 2),
    persen_bobot = round(persen_bobot, 2)
  )

persentase_ppd_bobot

# ============================================================
# PIE CHART PREVALENSI DEPRESI POSTPARTUM DENGAN PEMBOBOT
# ============================================================

grafik_ppd_bobot <- ggplot2::ggplot(
  persentase_ppd_bobot,
  ggplot2::aes(
    x = "",
    y = persen_bobot,
    fill = status
  )
) +
  ggplot2::geom_col(
    width = 1,
    color = "white"
  ) +
  ggplot2::coord_polar(
    theta = "y"
  ) +
  ggplot2::geom_text(
    ggplot2::aes(
      label = paste0(
        format(
          persen_bobot,
          decimal.mark = ",",
          nsmall = 2
        ),
        "%"
      )
    ),
    position = ggplot2::position_stack(vjust = 0.5),
    size = 5
  ) +
  ggplot2::labs(
    title = "Prevalensi Depresi Postpartum",
    fill = NULL
  ) +
  ggplot2::theme_void() +
  ggplot2::theme(
    legend.position = "bottom",
    plot.title = ggplot2::element_text(
      hjust = 0.5,
      face = "bold"
    )
  )

grafik_ppd_bobot

6.2 Visualisasi Variabel X

(bikin 1 chunk untuk mempermudah visualisasi variabel lain,

# ============================================================
# FUNGSI PROFILING LENGKAP PER VARIABEL X
# ============================================================

buat_profiling_x_lengkap <- function(data, x_var, nama_variabel, label_kategori) {
  
 data_plot <- data %>%
  dplyr::mutate(
    kode = as.character(haven::zap_labels(.data[[x_var]])),
    label = dplyr::recode(kode, !!!label_kategori),
    label = ifelse(is.na(label), "Missing", label),
    
    # Urutan kategori mengikuti urutan label_kategori
    label = factor(
      label,
      levels = c(unname(label_kategori), "Missing")
    ),
    
    status_label = dplyr::case_when(
      status_depresi == 0 ~ "Tidak PPD",
      status_depresi == 1 ~ "PPD",
      TRUE ~ "Missing"
    ),
    
    # Urutan status depresi juga dibuat tetap
    status_label = factor(
      status_label,
      levels = c("Tidak PPD", "PPD", "Missing")
    )
  )
  
  # Tabel kategori
  tabel_kategori <- data_plot %>%
    dplyr::count(kode, label, name = "n") %>%
    dplyr::mutate(
      nama_variabel = nama_variabel,
      variabel_x = x_var,
      persen = round(100 * n / sum(n), 2)
    ) %>%
    dplyr::select(nama_variabel, variabel_x, kode, label, n, persen)
  
  # Tabel bivariat
  tabel_bivariat <- data_plot %>%
    dplyr::count(kode, label, status_depresi, status_label, name = "n") %>%
    dplyr::group_by(kode, label) %>%
    dplyr::mutate(
      persen_dalam_kategori = round(100 * n / sum(n), 2)
    ) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(
      nama_variabel = nama_variabel,
      variabel_x = x_var
    ) %>%
    dplyr::select(
      nama_variabel,
      variabel_x,
      kode,
      label,
      status_depresi,
      status_label,
      n,
      persen_dalam_kategori
    )
  
  # Grafik distribusi kategori
  grafik_kategori <- ggplot2::ggplot(
    tabel_kategori,
    ggplot2::aes(x = label, y = persen)
  ) +
    ggplot2::geom_col() +
    ggplot2::geom_text(
      ggplot2::aes(label = paste0(n, " (", persen, "%)")),
      vjust = -0.3,
      size = 3.5
    ) +
    ggplot2::labs(
      title = paste("Distribusi", nama_variabel),
      x = nama_variabel,
      y = "Persentase"
    ) +
    ggplot2::theme_minimal() +
    ggplot2::theme(
      plot.title = ggplot2::element_text(hjust = 0.5, face = "bold")
    )
  
  # Grafik proporsi status depresi
  grafik_bivariat <- ggplot2::ggplot(
    tabel_bivariat,
    ggplot2::aes(
      x = label,
      y = persen_dalam_kategori / 100,
      fill = status_label
    )
  ) +
    ggplot2::geom_bar(stat = "identity") +
    ggplot2::geom_text(
      ggplot2::aes(
        label = scales::percent(persen_dalam_kategori / 100, accuracy = 0.01)
      ),
      position = ggplot2::position_stack(vjust = 0.5),
      size = 4
    ) +
    ggplot2::scale_y_continuous(labels = scales::percent) +
    ggplot2::labs(
      title = paste("Persentase Status Depresi Berdasarkan", nama_variabel),
      x = nama_variabel,
      y = "Persentase",
      fill = "Status depresi"
    ) +
    ggplot2::theme_minimal() +
    ggplot2::theme(
      plot.title = ggplot2::element_text(hjust = 0.5, face = "bold")
    )
  
  list(
    tabel_kategori = tabel_kategori,
    tabel_bivariat = tabel_bivariat,
    grafik_kategori = grafik_kategori,
    grafik_bivariat = grafik_bivariat
  )
}

6.2.1 X1_PAR

# ============================================================
# PROFILING X1: PARITAS 
# ============================================================

profiling_X1_PAR <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X1_PAR",
  nama_variabel = "Paritas",
  label_kategori = c(
    "0" = ">=2",
    "1" = "1"
  )
)

# Tampilkan tabel dan grafik
profiling_X1_PAR$tabel_kategori

profiling_X1_PAR$tabel_bivariat

profiling_X1_PAR$grafik_kategori

profiling_X1_PAR$grafik_bivariat

### X2_MP

# ============================================================
# PROFILING X2: METODE PERSALINAN
# ============================================================

# X2:
# 0 = Normal
# 1 = Selain normal
profiling_X2_MP <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X2_MP",
  nama_variabel = "Metode Persalinan",
  label_kategori = c(
    "0" = "Normal",
    "1" = "Selain normal"
  )
)

# Tampilkan tabel dan grafik
profiling_X2_MP$tabel_kategori

profiling_X2_MP$tabel_bivariat

profiling_X2_MP$grafik_kategori

profiling_X2_MP$grafik_bivariat

6.2.2 X3_KMN

# X3:
# 0 = Tidak ada komplikasi nifas
# 1 = Ada komplikasi nifas
profiling_X3_KMN <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X3_KMN",
  nama_variabel = "Komplikasi Masa Nifas",
  label_kategori = c(
    "0" = "Tidak Mengalami Komplikasi",
    "1" = "Mengalami Komplikasi"
  )
)

# Tampilkan tabel dan grafik
profiling_X3_KMN$tabel_kategori

profiling_X3_KMN$tabel_bivariat

profiling_X3_KMN$grafik_kategori

profiling_X3_KMN$grafik_bivariat

### X4_SKH

# X5:
# 0 = Diinginkan saat itu
# 1 = Diinginkan kemudian / tidak diinginkan
profiling_X4_SKH <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X4_SKH",
  nama_variabel = "Status Kehamilan",
  label_kategori = c(
    "0" = "Diinginkan saat itu",
    "1" = "Diinginkan kemudian / tidak diinginkan"
  )
)

profiling_X4_SKH$tabel_kategori

profiling_X4_SKH$tabel_bivariat

profiling_X4_SKH$grafik_kategori

profiling_X4_SKH$grafik_bivariat

X5_UI + X6_PI+ X7_WI + X8_JART + X9_JB,

6.2.3 X5_UI

# X6:
# 0 = Umur ibu 20-34 tahun
# 1 = Umur ibu <20 tahun atau >=35 tahun
profiling_X5_UI <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X5_UI",
  nama_variabel = "Umur Ibu",
  label_kategori = c(
    "0" = "20-34 tahun",
    "1" = "<20 tahun atau >=35 tahun"
  )
)

profiling_X5_UI$tabel_kategori

profiling_X5_UI$tabel_bivariat

profiling_X5_UI$grafik_kategori

profiling_X5_UI$grafik_bivariat

### X6_PI

# X8:
# 0 = SLTA atau lebih tinggi
# 1 = SMP atau lebih rendah
profiling_X6_PI <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X6_PI",
  nama_variabel = "Pendidikan Ibu",
  label_kategori = c(
    "0" = "SLTA atau lebih tinggi",
    "1" = "SMP atau lebih rendah"
  )
)

profiling_X6_PI$tabel_kategori

profiling_X6_PI$tabel_bivariat

profiling_X6_PI$grafik_kategori

profiling_X6_PI$grafik_bivariat

### X7_WI

# X9:
# 0 = Menengah-kaya (Q3-Q5)
# 1 = Miskin (Q1-Q2)
profiling_X7_WI <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X7_WI",
  nama_variabel = "Status Ekonomi",
  label_kategori = c(
    "0" = "Menengah-kaya (Q3-Q5)",
    "1" = "Miskin (Q1-Q2)"
  )
)
profiling_X7_WI$tabel_kategori

profiling_X7_WI$tabel_bivariat

profiling_X7_WI$grafik_kategori

profiling_X7_WI$grafik_bivariat

### X8_JART

# X10:
# 0 = 1-4 ART
# 1 = >=5 ART
profiling_X8_JART <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X8_JART",
  nama_variabel = "Jumlah Anggota Ruta",
  label_kategori = c(
    "0" = "1-4 ART",
    "1" = ">=5 ART"
  )
)

profiling_X8_JART$tabel_kategori

profiling_X8_JART$tabel_bivariat

profiling_X8_JART$grafik_kategori

profiling_X8_JART$grafik_bivariat

### X9_JB

# X11:
# 0 = 1 balita
# 1 = >=2 balita
profiling_X9_JB <- buat_profiling_x_lengkap(
  data = data_analisis,
  x_var = "X9_JB",
  nama_variabel = "Jumlah Balita dalam Ruta",
  label_kategori = c(
    "0" = "1 balita",
    "1" = ">=2 balita"
  )
)

profiling_X9_JB$tabel_kategori

profiling_X9_JB$tabel_bivariat

profiling_X9_JB$grafik_kategori

profiling_X9_JB$grafik_bivariat

7 Cek VIF

# ============================================================
# DATA UNTUK CEK VIF MODEL UTAMA
# ============================================================

library(dplyr)
library(haven)

data_vif <- data_analisis %>%
  dplyr::select(
    status_depresi,
    X1_PAR, X2_MP, X3_KMN, X4_SKH, X5_UI,
    X6_PI, X7_WI, X8_JART,
    X9_JB
  ) %>%
  dplyr::mutate(
    status_depresi = as.numeric(haven::zap_labels(status_depresi)),
    
    dplyr::across(
      .cols = -status_depresi,
      .fns = ~ factor(as.numeric(haven::zap_labels(.x)))
    )
  )

# ============================================================
# MODEL UNTUK CEK VIF
# ============================================================

model_vif <- glm(
  status_depresi ~ X1_PAR + X2_MP + X3_KMN + X4_SKH + X5_UI + X6_PI+ X7_WI + X8_JART + X9_JB,
  data = data_vif,
  family = binomial()
)

# ============================================================
# CEK VIF
# ============================================================

if (!require(car)) {
  install.packages("car")
  library(car)
}

## Loading required package: car

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:dplyr':
## 
##     recode

vif_model <- car::vif(model_vif)

vif_model

##   X1_PAR    X2_MP   X3_KMN   X4_SKH    X5_UI    X6_PI    X7_WI  X8_JART 
## 1.533971 1.074588 1.018061 1.186645 1.141252 1.112436 1.086592 1.271930 
##    X9_JB 
## 1.648822

Keterangan VIF < 5 = aman VIF 5–10 = mulai perlu diperhatikan VIF > 10 = indikasi multikolinearitas kuat

8 Odds Ratio Murni

library(haven)
library(dplyr)

# Lokasi data SPSS
data_model <- read.csv("C:/Users/nisri/Documents/KULIAH/SEMESTER 7/ZKRIPZWIT/POSTPARTUM DEPRESSION/OLAH DATA/After Sempro/SENSUS EKONOMI PROGRES/data_model4.csv")

# Membaca data
data_model 

# Variabel yang digunakan
variabel_x <- c(
  "X1", "X2", "X3", "X4", "X5",
  "X6", "X7", "X8", "X9"
)

# Fungsi menghitung crude odds ratio
hitung_or <- function(nama_x) {

  data_sementara <- data.frame(
    kategori = as_factor(data_model[[nama_x]]),
    status_depresi = as.numeric(data_model$status_depresi)
  ) %>%
    filter(
      !is.na(kategori),
      status_depresi %in% c(0, 1)
    )

  hasil <- data_sementara %>%
    group_by(kategori) %>%
    summarise(
      tidak_depresi = sum(status_depresi == 0),
      depresi = sum(status_depresi == 1),
      odds = depresi / tidak_depresi,
      .groups = "drop"
    ) %>%
    mutate(
      variabel = nama_x,
      odds_ratio = odds / first(odds),
      keterangan = ifelse(
        row_number() == 1,
        "Referensi",
        ""
      ),
      .before = 1
    )

  hasil
}

# Menghitung OR seluruh variabel
hasil_or <- bind_rows(
  lapply(variabel_x, hitung_or)
)

# Membulatkan hasil
hasil_or <- hasil_or %>%
  mutate(
    odds = round(odds, 3),
    odds_ratio = round(odds_ratio, 3)
  )

# Melihat hasil
hasil_or

Pemodelan

Buat semua var X bertipe Kategorik selain X5_NUM

# ============================================================
# MEMBUAT DATA MODEL: Y DAN X SIAP UNTUK REGRESI
# ============================================================

var_model_fix <- c(
  "status_depresi",
  "X1_PAR",
  "X2_MP",
  "X3_KMN",
  "X4_SKH",
  "X5_UI",
  "X5_NUM",
  "X6_PI",
  "X7_WI",
  "X8_JART",
  "X9_JB"
)

data_model_fix <- data_analisis %>%
  dplyr::select(dplyr::all_of(var_model_fix)) %>%
  dplyr::mutate(
    # Y tetap numeric 0/1 untuk regresi logistik
    status_depresi = as.numeric(haven::zap_labels(status_depresi)),
    
   # Semua X menjadi factor, KECUALI X5_num
    dplyr::across(
      .cols = -c(status_depresi, X5_NUM),
      .fns = ~ factor(as.numeric(haven::zap_labels(.x)))
    )
  )

# ============================================================
# CEK TIPE DATA SETELAH JADI DATA MODEL
# ============================================================

tipe_data_model_fix <- data.frame(
  variabel = names(data_model_fix),
  class = sapply(data_model_fix, function(x) paste(class(x), collapse = ", ")),
  typeof = sapply(data_model_fix, typeof),
  missing = sapply(data_model_fix, function(x) sum(is.na(x))),
  jumlah_kategori = sapply(data_model_fix, function(x) length(unique(na.omit(x))))
)

print(tipe_data_model_fix, row.names = FALSE)

##        variabel   class  typeof missing jumlah_kategori
##  status_depresi numeric  double       0               2
##          X1_PAR  factor integer       0               2
##           X2_MP  factor integer       0               2
##          X3_KMN  factor integer       0               2
##          X4_SKH  factor integer       0               2
##           X5_UI  factor integer       0               2
##          X5_NUM numeric  double       0              31
##           X6_PI  factor integer       0               2
##           X7_WI  factor integer       0               2
##         X8_JART  factor integer       0               2
##           X9_JB  factor integer       0               2

Seluruh Variabel

View(data_model_fix)

## Model 4 Var Signifikan of 9 (terbanyak signif)
model_fix <- logistf(
  formula = status_depresi ~ X1_PAR + X2_MP + X3_KMN + X4_SKH + X5_UI + X6_PI+ X7_WI + X8_JART + X9_JB,
  data = data_model_fix,
  pl = TRUE,
  firth = TRUE
)

summary(model_fix)

## logistf(formula = status_depresi ~ X1_PAR + X2_MP + X3_KMN + 
##     X4_SKH + X5_UI + X6_PI + X7_WI + X8_JART + X9_JB, data = data_model_fix, 
##     pl = TRUE, firth = TRUE)
## 
## Model fitted by Penalized ML
## Coefficients:
##                   coef  se(coef)  lower 0.95 upper 0.95      Chisq            p
## (Intercept) -4.9908590 0.5641617 -6.23862181 -3.9138579        Inf 0.0000000000
## X1_PAR1      1.2208482 0.5432018  0.06523094  2.3459600  4.2651229 0.0389024350
## X2_MP1      -0.2653435 0.4578746 -1.28542407  0.6336138  0.3129952 0.5758480178
## X3_KMN1      0.2128038 0.5276374 -1.00039255  1.2138865  0.1426026 0.7057072642
## X4_SKH1      1.0551830 0.4243110  0.16779105  1.9232328  5.3748475 0.0204291269
## X5_UI1      -0.4079788 0.5014098 -1.53474204  0.5716039  0.6266423 0.4285896365
## X6_PI1       0.9174327 0.3980782  0.10133781  1.7540587  4.8473528 0.0276885952
## X7_WI1       0.6932489 0.4794284 -0.35599787  1.6316360  1.7633076 0.1842124949
## X8_JART1     0.2576264 0.4244532 -0.62778102  1.1379738  0.3306250 0.5652913226
## X9_JB1       1.8661571 0.4854478  0.88966190  2.9325580 14.4498079 0.0001439445
##             method
## (Intercept)      2
## X1_PAR1          2
## X2_MP1           2
## X3_KMN1          2
## X4_SKH1          2
## X5_UI1           2
## X6_PI1           2
## X7_WI1           2
## X8_JART1         2
## X9_JB1           2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=35.05228 on 9 df, p=5.83314e-05, n=733
## Wald test = 220.8027 on 9 df, p = 0

exp(cbind(OR = coef(model_fix), confint(model_fix)))

##                      OR   Lower 95%   Upper 95%
## (Intercept) 0.006799821 0.001952545  0.01996334
## X1_PAR1     3.390061964 1.067405500 10.44329378
## X2_MP1      0.766942484 0.276533289  1.88440825
## X3_KMN1     1.237141877 0.367735059  3.36654318
## X4_SKH1     2.872500790 1.182689458  6.84304486
## X5_UI1      0.664992959 0.215511276  1.77110539
## X6_PI1      2.502856655 1.106650420  5.77800612
## X7_WI1      2.000203501 0.700474112  5.11223131
## X8_JART1    1.293855297 0.533774923  3.12043946
## X9_JB1      6.463410355 2.434306470 18.77559650

Hanya Variabel yang signifikan

model_2 <- logistf(
  formula = status_depresi ~ X1_PAR + X4_SKH + X6_PI + X9_JB,
  data = data_model_fix,
  pl = TRUE,
  firth = TRUE
)

summary(model_2)

## logistf(formula = status_depresi ~ X1_PAR + X4_SKH + X6_PI + 
##     X9_JB, data = data_model_fix, pl = TRUE, firth = TRUE)
## 
## Model fitted by Penalized ML
## Coefficients:
##                  coef  se(coef)  lower 0.95 upper 0.95     Chisq            p
## (Intercept) -5.029374 0.5118435 -6.14952715  -4.088200       Inf 0.000000e+00
## X1_PAR1      1.196753 0.5476590  0.06651979   2.295745  4.282113 3.851543e-02
## X4_SKH1      1.059640 0.4256919  0.18951656   1.905107  5.620792 1.774862e-02
## X6_PI1       1.012489 0.3986265  0.21881800   1.829267  6.229453 1.256426e-02
## X9_JB1       1.919830 0.4801997  0.98949202   2.937539 16.943981 3.849919e-05
##             method
## (Intercept)      2
## X1_PAR1          2
## X4_SKH1          2
## X6_PI1           2
## X9_JB1           2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=31.97775 on 4 df, p=1.933238e-06, n=733
## Wald test = 219.6243 on 4 df, p = 0

exp(cbind(OR = coef(model_2), confint(model_2)))

##                      OR   Lower 95%   Upper 95%
## (Intercept) 0.006542903 0.002134491  0.01676939
## X1_PAR1     3.309354633 1.068782116  9.93183661
## X4_SKH1     2.885332965 1.208665140  6.72012455
## X6_PI1      2.752444435 1.244604743  6.22931961
## X9_JB1      6.819796473 2.689867731 18.86935090

Seluruh Variabel (X5 Numerik)

model_3 <- logistf(
  formula = status_depresi ~X1_PAR + X2_MP + X3_KMN + X4_SKH +X5_NUM + X6_PI+ X7_WI + X8_JART + X9_JB,
  data = data_model_fix,
  pl = TRUE,
  firth = TRUE
)

summary(model_3)

## logistf(formula = status_depresi ~ X1_PAR + X2_MP + X3_KMN + 
##     X4_SKH + X5_NUM + X6_PI + X7_WI + X8_JART + X9_JB, data = data_model_fix, 
##     pl = TRUE, firth = TRUE)
## 
## Model fitted by Penalized ML
## Coefficients:
##                    coef  se(coef)   lower 0.95 upper 0.95       Chisq
## (Intercept) -3.00361835 1.4780653 -6.103691847 0.07329207  3.66102624
## X1_PAR1      0.70432647 0.6519980 -0.660822830 2.06570432  1.03921852
## X2_MP1      -0.17711776 0.4612183 -1.202675261 0.72977558  0.13592002
## X3_KMN1      0.16539984 0.5312323 -1.057383621 1.17199914  0.08562645
## X4_SKH1      0.99362778 0.4191595  0.116996178 1.84720703  4.89575258
## X5_NUM      -0.06452885 0.0438185 -0.158404742 0.02497390  1.97377454
## X6_PI1       0.81232649 0.3996926 -0.009276142 1.64819567  3.75595555
## X7_WI1       0.67154576 0.4815883 -0.381359382 1.61364281  1.64831308
## X8_JART1     0.37819537 0.4293099 -0.521355182 1.26661478  0.69145604
## X9_JB1       1.74888083 0.4906261  0.765965487 2.82446761 12.52389805
##                        p method
## (Intercept) 0.0556993155      2
## X1_PAR1     0.3080032935      2
## X2_MP1      0.7123712310      2
## X3_KMN1     0.7698126491      2
## X4_SKH1     0.0269228366      2
## X5_NUM      0.1600477852      2
## X6_PI1      0.0526197115      2
## X7_WI1      0.1991886852      2
## X8_JART1    0.4056695694      2
## X9_JB1      0.0004017798      2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=36.45338 on 9 df, p=3.293462e-05, n=733
## Wald test = 218.346 on 9 df, p = 0

exp(cbind(OR = coef(model_3), confint(model_3)))

##                     OR   Lower 95% Upper 95%
## (Intercept) 0.04960725 0.002234603  1.076045
## X1_PAR1     2.02248402 0.516426228  7.890854
## X2_MP1      0.83768113 0.300389516  2.074615
## X3_KMN1     1.17986478 0.347363457  3.228440
## X4_SKH1     2.70101541 1.124115134  6.342082
## X5_NUM      0.93750907 0.853504263  1.025288
## X6_PI1      2.25314381 0.990766748  5.197593
## X7_WI1      1.95726045 0.682932412  5.021069
## X8_JART1    1.45964809 0.593715410  3.548819
## X9_JB1      5.74816589 2.151070202 16.851971

9 Pengujian Kesesuaian Model

library(ResourceSelection)

## Warning: package 'ResourceSelection' was built under R version 4.3.3

## ResourceSelection 0.3-6   2023-06-27

# ============================================================
# UJI KESESUAIAN MODEL
# HOSMER-LEMESHOW GOODNESS OF FIT TEST
# ============================================================

library(ResourceSelection)

# Probabilitas prediksi dari model Firth
prob_prediksi <- predict(
  model_fix,
  type = "response"
)

# Hosmer-Lemeshow Test
uji_hl <- ResourceSelection::hoslem.test(
  x = model_fix$y,
  y = prob_prediksi,
  g = 10
)

uji_hl

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model_fix$y, prob_prediksi
## X-squared = 7.0188, df = 8, p-value = 0.5346

hasil_hosmer_lemeshow <- data.frame(
  Chi_Square = round(as.numeric(uji_hl$statistic), 3),
  df = as.numeric(uji_hl$parameter),
  p_value = round(uji_hl$p.value, 4)
)

hasil_hosmer_lemeshow

10 Uji Simultan

# ============================================================
# 10. UJI SIMULTAN
# PENALIZED LIKELIHOOD RATIO TEST
# ============================================================

uji_simultan <- logistf::logistftest(model_fix)

uji_simultan

## logistf::logistftest(object = model_fix)
## Model fitted by Penalized ML 
## 
## Factors fixed as follows:
## (Intercept)     X1_PAR1      X2_MP1     X3_KMN1     X4_SKH1      X5_UI1 
##          NA           0           0           0           0           0 
##      X6_PI1      X7_WI1    X8_JART1      X9_JB1 
##           0           0           0           0 
## 
## Likelihoods:
## Restricted model       Full model       difference 
##       -107.18385        -89.65771         17.52614 
## 
## Likelihood ratio test=35.05228 on 9 df, p=5.83314e-05

11 Uji Parsial

11.1 Seluruh Variabel

View(data_model_fix)

## Model 4 Var Signifikan of 9 (terbanyak signif)
model_fix <- logistf(
  formula = status_depresi ~ X1_PAR + X2_MP + X3_KMN + X4_SKH + X5_UI + X6_PI+ X7_WI + X8_JART + X9_JB,
  data = data_model_fix,
  pl = TRUE,
  firth = TRUE
)

summary(model_fix)

## logistf(formula = status_depresi ~ X1_PAR + X2_MP + X3_KMN + 
##     X4_SKH + X5_UI + X6_PI + X7_WI + X8_JART + X9_JB, data = data_model_fix, 
##     pl = TRUE, firth = TRUE)
## 
## Model fitted by Penalized ML
## Coefficients:
##                   coef  se(coef)  lower 0.95 upper 0.95      Chisq            p
## (Intercept) -4.9908590 0.5641617 -6.23862181 -3.9138579        Inf 0.0000000000
## X1_PAR1      1.2208482 0.5432018  0.06523094  2.3459600  4.2651229 0.0389024350
## X2_MP1      -0.2653435 0.4578746 -1.28542407  0.6336138  0.3129952 0.5758480178
## X3_KMN1      0.2128038 0.5276374 -1.00039255  1.2138865  0.1426026 0.7057072642
## X4_SKH1      1.0551830 0.4243110  0.16779105  1.9232328  5.3748475 0.0204291269
## X5_UI1      -0.4079788 0.5014098 -1.53474204  0.5716039  0.6266423 0.4285896365
## X6_PI1       0.9174327 0.3980782  0.10133781  1.7540587  4.8473528 0.0276885952
## X7_WI1       0.6932489 0.4794284 -0.35599787  1.6316360  1.7633076 0.1842124949
## X8_JART1     0.2576264 0.4244532 -0.62778102  1.1379738  0.3306250 0.5652913226
## X9_JB1       1.8661571 0.4854478  0.88966190  2.9325580 14.4498079 0.0001439445
##             method
## (Intercept)      2
## X1_PAR1          2
## X2_MP1           2
## X3_KMN1          2
## X4_SKH1          2
## X5_UI1           2
## X6_PI1           2
## X7_WI1           2
## X8_JART1         2
## X9_JB1           2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=35.05228 on 9 df, p=5.83314e-05, n=733
## Wald test = 220.8027 on 9 df, p = 0

exp(cbind(OR = coef(model_fix), confint(model_fix)))

##                      OR   Lower 95%   Upper 95%
## (Intercept) 0.006799821 0.001952545  0.01996334
## X1_PAR1     3.390061964 1.067405500 10.44329378
## X2_MP1      0.766942484 0.276533289  1.88440825
## X3_KMN1     1.237141877 0.367735059  3.36654318
## X4_SKH1     2.872500790 1.182689458  6.84304486
## X5_UI1      0.664992959 0.215511276  1.77110539
## X6_PI1      2.502856655 1.106650420  5.77800612
## X7_WI1      2.000203501 0.700474112  5.11223131
## X8_JART1    1.293855297 0.533774923  3.12043946
## X9_JB1      6.463410355 2.434306470 18.77559650

11.2 Hanya Variabel yang signifikan

model_2 <- logistf(
  formula = status_depresi ~ X1_PAR + X4_SKH + X6_PI + X9_JB,
  data = data_model_fix,
  pl = TRUE,
  firth = TRUE
)

summary(model_2)

## logistf(formula = status_depresi ~ X1_PAR + X4_SKH + X6_PI + 
##     X9_JB, data = data_model_fix, pl = TRUE, firth = TRUE)
## 
## Model fitted by Penalized ML
## Coefficients:
##                  coef  se(coef)  lower 0.95 upper 0.95     Chisq            p
## (Intercept) -5.029374 0.5118435 -6.14952715  -4.088200       Inf 0.000000e+00
## X1_PAR1      1.196753 0.5476590  0.06651979   2.295745  4.282113 3.851543e-02
## X4_SKH1      1.059640 0.4256919  0.18951656   1.905107  5.620792 1.774862e-02
## X6_PI1       1.012489 0.3986265  0.21881800   1.829267  6.229453 1.256426e-02
## X9_JB1       1.919830 0.4801997  0.98949202   2.937539 16.943981 3.849919e-05
##             method
## (Intercept)      2
## X1_PAR1          2
## X4_SKH1          2
## X6_PI1           2
## X9_JB1           2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=31.97775 on 4 df, p=1.933238e-06, n=733
## Wald test = 219.6243 on 4 df, p = 0

exp(cbind(OR = coef(model_2), confint(model_2)))

##                      OR   Lower 95%   Upper 95%
## (Intercept) 0.006542903 0.002134491  0.01676939
## X1_PAR1     3.309354633 1.068782116  9.93183661
## X4_SKH1     2.885332965 1.208665140  6.72012455
## X6_PI1      2.752444435 1.244604743  6.22931961
## X9_JB1      6.819796473 2.689867731 18.86935090

11.3 Seluruh Variabel (X5 Numerik)

model_3 <- logistf(
  formula = status_depresi ~X1_PAR + X2_MP + X3_KMN + X4_SKH +X5_NUM + X6_PI+ X7_WI + X8_JART + X9_JB,
  data = data_model_fix,
  pl = TRUE,
  firth = TRUE
)

summary(model_3)

## logistf(formula = status_depresi ~ X1_PAR + X2_MP + X3_KMN + 
##     X4_SKH + X5_NUM + X6_PI + X7_WI + X8_JART + X9_JB, data = data_model_fix, 
##     pl = TRUE, firth = TRUE)
## 
## Model fitted by Penalized ML
## Coefficients:
##                    coef  se(coef)   lower 0.95 upper 0.95       Chisq
## (Intercept) -3.00361835 1.4780653 -6.103691847 0.07329207  3.66102624
## X1_PAR1      0.70432647 0.6519980 -0.660822830 2.06570432  1.03921852
## X2_MP1      -0.17711776 0.4612183 -1.202675261 0.72977558  0.13592002
## X3_KMN1      0.16539984 0.5312323 -1.057383621 1.17199914  0.08562645
## X4_SKH1      0.99362778 0.4191595  0.116996178 1.84720703  4.89575258
## X5_NUM      -0.06452885 0.0438185 -0.158404742 0.02497390  1.97377454
## X6_PI1       0.81232649 0.3996926 -0.009276142 1.64819567  3.75595555
## X7_WI1       0.67154576 0.4815883 -0.381359382 1.61364281  1.64831308
## X8_JART1     0.37819537 0.4293099 -0.521355182 1.26661478  0.69145604
## X9_JB1       1.74888083 0.4906261  0.765965487 2.82446761 12.52389805
##                        p method
## (Intercept) 0.0556993155      2
## X1_PAR1     0.3080032935      2
## X2_MP1      0.7123712310      2
## X3_KMN1     0.7698126491      2
## X4_SKH1     0.0269228366      2
## X5_NUM      0.1600477852      2
## X6_PI1      0.0526197115      2
## X7_WI1      0.1991886852      2
## X8_JART1    0.4056695694      2
## X9_JB1      0.0004017798      2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=36.45338 on 9 df, p=3.293462e-05, n=733
## Wald test = 218.346 on 9 df, p = 0

exp(cbind(OR = coef(model_3), confint(model_3)))

##                     OR   Lower 95% Upper 95%
## (Intercept) 0.04960725 0.002234603  1.076045
## X1_PAR1     2.02248402 0.516426228  7.890854
## X2_MP1      0.83768113 0.300389516  2.074615
## X3_KMN1     1.17986478 0.347363457  3.228440
## X4_SKH1     2.70101541 1.124115134  6.342082
## X5_NUM      0.93750907 0.853504263  1.025288
## X6_PI1      2.25314381 0.990766748  5.197593
## X7_WI1      1.95726045 0.682932412  5.021069
## X8_JART1    1.45964809 0.593715410  3.548819
## X9_JB1      5.74816589 2.151070202 16.851971