1 1. Persiapan

Penting: Instalasi package dilakukan satu kali melalui Console RStudio, bukan saat Knit. Jalankan sekali di Console bila package belum tersedia:

install.packages(
  c("readxl", "dplyr", "tidyr", "ggplot2", "knitr", "rmarkdown", "survival"),
  repos = "https://cloud.r-project.org"
)

Setelah instalasi selesai, gunakan tombol Knit untuk membuat HTML.

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(knitr)
library(survival)

2 2. Import data

Letakkan file Data Penelitian Diolah.xlsx dalam folder kerja RStudio yang sama dengan file Rmd ini. Jika lokasinya berbeda, ubah nilai file_data.

file_data <- "Data Penelitian Diolah.xlsx"

if (!file.exists(file_data)) {
  stop("File Excel tidak ditemukan. Letakkan file Excel di folder kerja atau ubah nilai file_data.")
}

data_awal <- read_excel(file_data, sheet = "OLAH")
names(data_awal)
##  [1] "TAHUN"     "URUTAN"    "WEIGHT"    "KODE_PROV" "KODE_KAB"  "Y"        
##  [7] "X1"        "X2"        "X3"        "X4"        "X5"        "X6"       
## [13] "X7"
dim(data_awal)
## [1] 2135   13

3 3. Pemeriksaan struktur data

str(data_awal)
## tibble [2,135 × 13] (S3: tbl_df/tbl/data.frame)
##  $ TAHUN    : num [1:2135] 20258 20258 20258 20258 20258 ...
##  $ URUTAN   : num [1:2135] 546476 546713 547306 547310 547347 ...
##  $ WEIGHT   : num [1:2135] 45 93 204 187 69 127 100 153 129 89 ...
##  $ KODE_PROV: num [1:2135] 53 53 53 53 53 53 53 53 53 53 ...
##  $ KODE_KAB : num [1:2135] 21 2 19 19 21 2 2 21 6 18 ...
##  $ Y        : num [1:2135] 1 3 1 1 4 12 12 1 3 1 ...
##  $ X1       : num [1:2135] 1 1 0 1 0 0 0 1 0 1 ...
##  $ X2       : num [1:2135] 17 42 51 43 35 26 24 50 60 80 ...
##  $ X3       : num [1:2135] 0 1 1 1 1 0 0 1 1 1 ...
##  $ X4       : num [1:2135] 0 0 0 0 0 0 0 0 0 0 ...
##  $ X5       : num [1:2135] 0 1 1 1 1 0 0 1 0 0 ...
##  $ X6       : num [1:2135] 0 0 0 0 0 1 1 0 0 0 ...
##  $ X7       : num [1:2135] 0 2 0 0 0 1 1 0 0 0 ...
knitr::kable(
  data.frame(
    Variabel = names(data_awal),
    Missing = sapply(data_awal, function(x) sum(is.na(x))),
    Jumlah_nilai_unik = sapply(data_awal, function(x) dplyr::n_distinct(x))
  ),
  caption = "Pemeriksaan awal variabel"
)
Pemeriksaan awal variabel
Variabel Missing Jumlah_nilai_unik
TAHUN TAHUN 0 1
URUTAN URUTAN 0 2135
WEIGHT WEIGHT 0 361
KODE_PROV KODE_PROV 0 1
KODE_KAB KODE_KAB 0 22
Y Y 0 31
X1 X1 0 2
X2 X2 0 73
X3 X3 0 2
X4 X4 0 2
X5 X5 0 2
X6 X6 0 2
X7 X7 2 4

Dataset yang diperiksa menggunakan 2.135 observasi dan memuat variabel Y, X1–X7, serta WEIGHT. Kode berikut tetap memeriksa ulang kondisi file setiap kali dijalankan.

4 4. Menyiapkan variabel penelitian

Definisi kerja yang digunakan: - Y: waktu tunggu memperoleh pekerjaan, dalam bulan. - X1: jenis kelamin. - X2: umur, dalam tahun. - X3: status perkawinan. - X4: klasifikasi wilayah tempat tinggal. - X5: pengalaman kerja. - X6: pelatihan/kursus. - X7: pendidikan. - WEIGHT: bobot individu.

Periksa kembali buku kode Sakernas/kode pengolahan Anda untuk memastikan arti kategori 0, 1, dan kategori X7. Model memperlakukan X1, X3, X4, X5, X6, dan X7 sebagai variabel kategorik; X2 sebagai numerik.

variabel_wajib <- c("Y", "X1", "X2", "X3", "X4", "X5", "X6", "X7", "WEIGHT")
tidak_ada <- setdiff(variabel_wajib, names(data_awal))

if (length(tidak_ada) > 0) {
  stop(paste("Kolom berikut tidak ditemukan:", paste(tidak_ada, collapse = ", ")))
}

data_analisis <- data_awal %>%
  mutate(
    Y = as.numeric(Y),
    X2 = as.numeric(X2),
    WEIGHT = as.numeric(WEIGHT),
    across(c(X1, X3, X4, X5, X6, X7), as.factor)
  ) %>%
  filter(
    !is.na(Y), Y > 0,
    !is.na(WEIGHT), WEIGHT > 0
  )

cat("Jumlah observasi awal:", nrow(data_awal), "\n")
## Jumlah observasi awal: 2135
cat("Jumlah observasi setelah validasi Y dan WEIGHT:", nrow(data_analisis), "\n")
## Jumlah observasi setelah validasi Y dan WEIGHT: 2135
cat("Missing X7:", sum(is.na(data_analisis$X7)), "\n")
## Missing X7: 2

5 5. Pembersihan data untuk model

Model lengkap menggunakan observasi yang memiliki data lengkap pada seluruh variabel model. Observasi dengan X7 kosong tidak diisi dengan kategori buatan; observasi tersebut dikeluarkan dari model lengkap dan jumlahnya dilaporkan.

n_sebelum <- nrow(data_analisis)

data_model <- data_analisis %>%
  select(Y, X1, X2, X3, X4, X5, X6, X7, WEIGHT) %>%
  drop_na()

cat("Observasi sebelum menghapus missing kovariat:", n_sebelum, "\n")
## Observasi sebelum menghapus missing kovariat: 2135
cat("Observasi yang digunakan dalam model:", nrow(data_model), "\n")
## Observasi yang digunakan dalam model: 2133
cat("Observasi yang dikeluarkan karena missing:", n_sebelum - nrow(data_model), "\n")
## Observasi yang dikeluarkan karena missing: 2
if (nrow(data_model) < 30) {
  stop("Jumlah observasi lengkap terlalu sedikit untuk melanjutkan model.")
}

6 6. Statistik deskriptif

ringkasan_Y <- data_model %>%
  summarise(
    N = n(),
    Minimum = min(Y),
    Q1 = quantile(Y, 0.25),
    Median = median(Y),
    Mean = mean(Y),
    Q3 = quantile(Y, 0.75),
    Maksimum = max(Y),
    SD = sd(Y)
  )

knitr::kable(ringkasan_Y, digits = 2,
             caption = "Statistik deskriptif waktu tunggu (bulan)")
Statistik deskriptif waktu tunggu (bulan)
N Minimum Q1 Median Mean Q3 Maksimum SD
2133 1 1 2 4.23 5 60 5.96 
frekuensi_Y <- data_model %>%
  count(Y, name = "Frekuensi") %>%
  arrange(Y)

knitr::kable(frekuensi_Y, caption = "Distribusi frekuensi Y")
Distribusi frekuensi Y
Y Frekuensi
1 955
2 286
3 174
4 83
5 102
6 119
7 38
8 62
9 28
10 45
11 47
12 121
13 6
14 5
15 3
17 2
18 2
19 2
20 1
21 2
22 1
23 1
24 28
25 1
27 2
28 1
36 4
39 1
40 1
48 3
60 7 
ggplot(data_model, aes(x = Y)) +
  geom_histogram(binwidth = 1, boundary = 0.5, color = "white") +
  labs(
    title = "Distribusi Waktu Tunggu Memperoleh Pekerjaan",
    x = "Waktu tunggu (bulan)",
    y = "Frekuensi"
  ) +
  theme_minimal()

7 7. Distribusi variabel bebas

variabel_kategorik <- c("X1", "X3", "X4", "X5", "X6", "X7")

for (v in variabel_kategorik) {
  cat("\n\n### Distribusi", v, "\n\n")
  print(knitr::kable(as.data.frame(table(data_model[[v]])),
                     col.names = c("Kategori_kode", "Frekuensi")))
}
## 
## 
## ### Distribusi X1 
## 
## 
## 
## |Kategori_kode | Frekuensi|
## |:-------------|---------:|
## |0             |      1005|
## |1             |      1128|
## 
## 
## ### Distribusi X3 
## 
## 
## 
## |Kategori_kode | Frekuensi|
## |:-------------|---------:|
## |0             |       981|
## |1             |      1152|
## 
## 
## ### Distribusi X4 
## 
## 
## 
## |Kategori_kode | Frekuensi|
## |:-------------|---------:|
## |0             |      1607|
## |1             |       526|
## 
## 
## ### Distribusi X5 
## 
## 
## 
## |Kategori_kode | Frekuensi|
## |:-------------|---------:|
## |0             |      1328|
## |1             |       805|
## 
## 
## ### Distribusi X6 
## 
## 
## 
## |Kategori_kode | Frekuensi|
## |:-------------|---------:|
## |0             |      1774|
## |1             |       359|
## 
## 
## ### Distribusi X7 
## 
## 
## 
## |Kategori_kode | Frekuensi|
## |:-------------|---------:|
## |0             |      1026|
## |1             |       387|
## |2             |       720|
ggplot(data_model, aes(x = X2)) +
  geom_histogram(binwidth = 5, color = "white") +
  labs(title = "Distribusi Umur Responden",
       x = "Umur (tahun)", y = "Frekuensi") +
  theme_minimal()

8 8. Metode dan batas interpretasi

Y adalah durasi positif dalam bulan dan pada file ini tidak tersedia indikator sensor/censoring. Dengan demikian, model di bawah merupakan model AFT parametrik untuk durasi yang teramati/selesai, bukan analisis hazard kohort pencari kerja dengan observasi tersensor.

Karena responden berasal dari kelompok yang memenuhi syarat pertanyaan Sakernas dan sudah memperoleh pekerjaan, hasil harus ditafsirkan sebagai hubungan karakteristik dengan lama durasi pada kelompok responden yang teramati tersebut. Hasil tidak otomatis mewakili seluruh pencari kerja NTT atau peluang memperoleh pekerjaan bagi semua penganggur.

Model AFT lognormal dan Weibull diestimasi sebagai kandidat. Pemilihan kandidat menggunakan AIC pada sampel dan representasi likelihood yang sama. Bobot WEIGHT dimasukkan sebagai bobot kasus. Karena file olahan ini tidak memuat informasi strata dan PSU, standard error yang dihasilkan tidak memperhitungkan keseluruhan desain sampel kompleks Sakernas.

9 9. Estimasi model AFT

formula_aft <- survival::Surv(Y, rep(1, nrow(data_model))) ~
  X1 + X2 + X3 + X4 + X5 + X6 + X7

model_lognormal <- survival::survreg(
  formula_aft,
  data = data_model,
  dist = "lognormal",
  weights = WEIGHT
)

model_weibull <- survival::survreg(
  formula_aft,
  data = data_model,
  dist = "weibull",
  weights = WEIGHT
)

summary(model_lognormal)
## 
## Call:
## survival::survreg(formula = formula_aft, data = data_model, weights = WEIGHT, 
##     dist = "lognormal")
##                Value Std. Error      z       p
## (Intercept)  0.58027    0.00640  90.67 < 2e-16
## X11          0.03221    0.00363   8.87 < 2e-16
## X2           0.00106    0.00017   6.26 3.9e-10
## X31         -0.01434    0.00473  -3.03  0.0024
## X41          0.27940    0.00395  70.79 < 2e-16
## X51         -0.03354    0.00388  -8.65 < 2e-16
## X61         -0.06604    0.00504 -13.11 < 2e-16
## X71          0.71583    0.00545 131.23 < 2e-16
## X72          0.33043    0.00445  74.28 < 2e-16
## Log(scale)  -0.04353    0.00131 -33.31 < 2e-16
## 
## Scale= 0.957 
## 
## Log Normal distribution
## Loglik(model)= -687064.4   Loglik(intercept only)= -702822
##  Chisq= 31515.35 on 8 degrees of freedom, p= 0 
## Number of Newton-Raphson Iterations: 3 
## n= 2133
summary(model_weibull)
## 
## Call:
## survival::survreg(formula = formula_aft, data = data_model, weights = WEIGHT, 
##     dist = "weibull")
##                 Value Std. Error      z       p
## (Intercept)  1.086810   0.007303 148.81 < 2e-16
## X11          0.031932   0.003903   8.18 2.8e-16
## X2          -0.000263   0.000193  -1.36  0.1725
## X31          0.016422   0.005018   3.27  0.0011
## X41          0.343902   0.004277  80.40 < 2e-16
## X51         -0.086934   0.004164 -20.88 < 2e-16
## X61         -0.055896   0.005290 -10.57 < 2e-16
## X71          0.748109   0.005841 128.09 < 2e-16
## X72          0.365904   0.004795  76.31 < 2e-16
## Log(scale)   0.018843   0.001310  14.39 < 2e-16
## 
## Scale= 1.02 
## 
## Weibull distribution
## Loglik(model)= -726502.2   Loglik(intercept only)= -744439.2
##  Chisq= 35873.97 on 8 degrees of freedom, p= 0 
## Number of Newton-Raphson Iterations: 5 
## n= 2133

10 10. Perbandingan kandidat model

tabel_aic <- data.frame(
  Model = c("AFT Lognormal", "AFT Weibull"),
  AIC = c(AIC(model_lognormal), AIC(model_weibull))
) %>%
  arrange(AIC)

knitr::kable(tabel_aic, digits = 2,
             caption = "Perbandingan AIC model parametrik")
Perbandingan AIC model parametrik
Model AIC
AFT Lognormal 1374149
AFT Weibull 1453024 
model_terpilih <- if (AIC(model_lognormal) <= AIC(model_weibull)) {
  model_lognormal
} else {
  model_weibull
}

nama_model_terpilih <- if (identical(model_terpilih, model_lognormal)) {
  "AFT Lognormal"
} else {
  "AFT Weibull"
}

cat("Model dengan AIC lebih kecil:", nama_model_terpilih)
## Model dengan AIC lebih kecil: AFT Lognormal

11 11. Rasio waktu (time ratio)

Pada model AFT, koefisien positif menunjukkan waktu tunggu cenderung lebih panjang, sedangkan koefisien negatif menunjukkan waktu tunggu cenderung lebih pendek, dengan variabel lain konstan. Eksponensial koefisien merupakan time ratio.

koef <- summary(model_terpilih)$table
koef <- as.data.frame(koef)
koef$Variabel <- rownames(koef)
rownames(koef) <- NULL

nama_kolom <- names(koef)
kolom_est <- grep("Value", nama_kolom, value = TRUE)[1]
kolom_se <- grep("Std", nama_kolom, value = TRUE)[1]
kolom_p <- grep("p", nama_kolom, value = TRUE, ignore.case = TRUE)
kolom_p <- kolom_p[length(kolom_p)]

hasil_model <- koef %>%
  filter(Variabel != "Log(scale)") %>%
  mutate(
    Time_Ratio = exp(.data[[kolom_est]]),
    CI_bawah = exp(.data[[kolom_est]] - 1.96 * .data[[kolom_se]]),
    CI_atas = exp(.data[[kolom_est]] + 1.96 * .data[[kolom_se]])
  )

knitr::kable(
  hasil_model %>% select(Variabel, all_of(kolom_est), all_of(kolom_se),
                         Time_Ratio, CI_bawah, CI_atas),
  digits = 4,
  caption = "Estimasi koefisien dan time ratio"
)
Estimasi koefisien dan time ratio
Variabel Value Std. Error Time_Ratio CI_bawah CI_atas
(Intercept) 0.5803 0.0064 1.7865 1.7642 1.8091
X11 0.0322 0.0036 1.0327 1.0254 1.0401
X2 0.0011 0.0002 1.0011 1.0007 1.0014
X31 -0.0143 0.0047 0.9858 0.9767 0.9949
X41 0.2794 0.0039 1.3223 1.3121 1.3326
X51 -0.0335 0.0039 0.9670 0.9597 0.9744
X61 -0.0660 0.0050 0.9361 0.9269 0.9454
X71 0.7158 0.0055 2.0459 2.0241 2.0679
X72 0.3304 0.0044 1.3916 1.3795 1.4038

12 12. Interpretasi otomatis secara hati-hati

cat("Model terpilih berdasarkan AIC:", nama_model_terpilih, "\n")
## Model terpilih berdasarkan AIC: AFT Lognormal
cat("Jumlah observasi model:", nrow(data_model), "\n")
## Jumlah observasi model: 2133
cat("Median waktu tunggu:", median(data_model$Y), "bulan\n")
## Median waktu tunggu: 2 bulan
cat("Rata-rata waktu tunggu:", round(mean(data_model$Y), 2), "bulan\n")
## Rata-rata waktu tunggu: 4.23 bulan

Interpretasikan time ratio sebagai berikut: - Time ratio lebih besar dari 1: durasi waktu tunggu diperkirakan lebih panjang dibandingkan kategori referensi (untuk variabel kategorik), atau meningkat seiring kenaikan satu unit X2. - Time ratio lebih kecil dari 1: durasi waktu tunggu diperkirakan lebih pendek. - Time ratio sama dengan 1: tidak ada perbedaan durasi berdasarkan estimasi titik. - Signifikansi statistik perlu dilihat dari p-value dan interval kepercayaan, bukan dari time ratio saja.

13 13. Catatan akhir

Analisis ini bersifat kondisional pada responden yang masuk ke dalam sampel olahan dan memiliki Y positif serta kovariat lengkap. Jangan menyimpulkan pengaruh terhadap seluruh pencari kerja/penganggur di NTT tanpa data durasi yang mencakup individu yang belum memperoleh pekerjaan dan informasi censoring yang benar.

Bobot WEIGHT telah digunakan sebagai bobot kasus, tetapi tanpa variabel strata dan PSU, inferensi belum sepenuhnya mengakomodasi desain survei Sakernas. Untuk hasil penelitian final/publikasi, verifikasi kembali definisi populasi sasaran, aturan skip pertanyaan Y, kode kategori X1–X7, serta rancangan sampel dari dokumentasi Sakernas.