1 Pendahuluan

Regresi probit adalah salah satu metode yang digunakan untuk menggambarkan hubungan variabel respon yang bersifat biner dengan satu atau lebih variabel bebas yang bersifat kontinu, kategorik atau kombinasi keduanya (Agresti, 2002).

Model regresi probit merupakan pengembangan dari model regresi logistik dengan menggunakan fungsi normal kumulatif (Cummulatif Density Function Normal/ CDF Normal) sedangkan pada regresi logistik menggunakan fungsi logistik kumulatif. Istilah probit berasal dari singkatan probability unit yang dikenalkan pada tahun 1930-an oleh Chester Bliss.

Model regresi probit didefinisikan sebagai:

Ringkasnya adalah:

\[ P(Y = 1 | X) = \Phi(X\beta) \]

dimana:

  • \(P(Y = 1 | X)\): Probabilitas respon.

  • \(\Phi\): Fungsi distribusi kumulatif normal standar.

  • \(X\beta\): Kombinasi linier variabel independen dan parameter.

Pendugaan parameter diperoleh dengan menggunakan metode estimasi Maximum Likelihood (MLE), yang merepresentasikan probabilitas model yang diestimasi menghasilkan data yang diamati. Pengamatan diasumsikan saling bebas (independen).

library(readxl)
library(broom)
library(caret)
library(DataExplorer)
library(grid)
library(ISLR)
library(pscl)
library(tidyverse)
library(dplyr)

Load Data

glimpse(dataku)
## Rows: 5,551
## Columns: 17
## $ R101      <int> 76, 31, 81, 72, 32, 52, 16, 65, 52, 36, 73, 31, 72, 34, 62, …
## $ R102      <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ R105      <int> 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 1, …
## $ R401      <int> 5, 6, 4, 4, 4, 3, 4, 7, 4, 8, 5, 5, 5, 6, 3, 4, 7, 4, 5, 4, …
## $ R403      <int> 3, 6, 3, 3, 3, 3, 6, 6, 3, 6, 3, 6, 6, 6, 3, 3, 6, 3, 9, 3, …
## $ R404      <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ R405      <int> 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 2, 2, 1, 1, 2, 1, 2, 2, 1, 2, …
## $ R407      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ASI_Eks   <int> 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, …
## $ R1401     <int> 1, 4, 0, 5, 3, 0, 2, 5, 3, 2, 2, 5, 5, 1, 1, 5, 0, 0, 5, 1, …
## $ TAMAT_KRT <int> 4, 1, 2, 3, 3, 4, 4, 1, 2, 2, 3, 4, 2, 3, 2, 2, 2, 4, 1, 4, …
## $ KERJA_KRT <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ NKAPITA   <int> 3, 2, 4, 2, 1, 5, 2, 2, 4, 1, 1, 4, 1, 2, 5, 2, 4, 3, 1, 5, …
## $ MISKIN    <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ SEX_KRT   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, …
## $ HH_FAMILY <int> 2, 3, 2, 2, 2, 1, 2, 3, 2, 3, 2, 2, 2, 3, 1, 3, 3, 3, 2, 2, …
## $ FWT       <chr> "125.55", "34.02", "102.82", "91.84", "346.76", "106.3", "19…

Eksplore Data

dataku$ASI_Eks <- as.factor(dataku$ASI_Eks)
dataku$R105 <- as.factor(dataku$R105)
dataku$R405 <- as.factor(dataku$R405)
dataku$TAMAT_KRT <- as.factor(dataku$TAMAT_KRT)
dataku$KERJA_KRT <- as.factor(dataku$KERJA_KRT)
dataku$NKAPITA <- as.factor(dataku$NKAPITA)
dataku$SEX_KRT <- as.factor(dataku$SEX_KRT)
dataku$HH_FAMILY <- as.factor(dataku$HH_FAMILY)
dataku$MISKIN <- as.factor(dataku$MISKIN)
str(dataku)
## 'data.frame':    5551 obs. of  17 variables:
##  $ R101     : int  76 31 81 72 32 52 16 65 52 36 ...
##  $ R102     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ R105     : Factor w/ 2 levels "1","2": 1 1 1 2 2 2 2 2 2 2 ...
##  $ R401     : int  5 6 4 4 4 3 4 7 4 8 ...
##  $ R403     : int  3 6 3 3 3 3 6 6 3 6 ...
##  $ R404     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ R405     : Factor w/ 2 levels "1","2": 1 1 2 2 1 1 1 1 2 1 ...
##  $ R407     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ASI_Eks  : Factor w/ 2 levels "0","1": 2 2 2 1 2 2 2 1 2 1 ...
##  $ R1401    : int  1 4 0 5 3 0 2 5 3 2 ...
##  $ TAMAT_KRT: Factor w/ 4 levels "1","2","3","4": 4 1 2 3 3 4 4 1 2 2 ...
##  $ KERJA_KRT: Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
##  $ NKAPITA  : Factor w/ 5 levels "1","2","3","4",..: 3 2 4 2 1 5 2 2 4 1 ...
##  $ MISKIN   : Factor w/ 2 levels "1","2": 2 2 2 2 2 2 2 2 2 1 ...
##  $ SEX_KRT  : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ HH_FAMILY: Factor w/ 3 levels "1","2","3": 2 3 2 2 2 1 2 3 2 3 ...
##  $ FWT      : chr  "125.55" "34.02" "102.82" "91.84" ...

Tampak beberapa variabel yang dianggap berpotensi mempengaruhi ASI Eksklusif, sekarang sudah bertipe factor.

summary(dataku)
##       R101            R102       R105          R401             R403     
##  Min.   :11.00   Min.   : 1.00   1:2107   Min.   : 2.000   Min.   :3.00  
##  1st Qu.:18.00   1st Qu.: 4.00   2:3444   1st Qu.: 4.000   1st Qu.:3.00  
##  Median :35.00   Median : 9.00            Median : 5.000   Median :3.00  
##  Mean   :44.36   Mean   :19.85            Mean   : 5.133   Mean   :3.94  
##  3rd Qu.:71.00   3rd Qu.:20.00            3rd Qu.: 6.000   3rd Qu.:6.00  
##  Max.   :97.00   Max.   :79.00            Max.   :19.000   Max.   :9.00  
##       R404   R405          R407   ASI_Eks      R1401       TAMAT_KRT KERJA_KRT
##  Min.   :1   1:2898   Min.   :0   0:1601   Min.   :0.000   1: 662    0: 274   
##  1st Qu.:1   2:2653   1st Qu.:0   1:3950   1st Qu.:1.000   2:1502    1:5277   
##  Median :1            Median :0            Median :3.000   3:1044             
##  Mean   :1            Mean   :0            Mean   :2.682   4:2343             
##  3rd Qu.:1            3rd Qu.:0            3rd Qu.:4.000                      
##  Max.   :1            Max.   :0            Max.   :5.000                      
##  NKAPITA  MISKIN   SEX_KRT  HH_FAMILY     FWT           
##  1:1629   1: 895   1:5196   1: 349    Length:5551       
##  2:1308   2:4656   2: 355   2:3100    Class :character  
##  3:1053                     3:2102    Mode  :character  
##  4: 895                                                 
##  5: 666                                                 
##

Pemeriksaan distribusi variabel-variabel. 1

plot_bar(dataku[,c("ASI_Eks","R105","R405","TAMAT_KRT","KERJA_KRT","NKAPITA", 
                      "SEX_KRT","HH_FAMILY","MISKIN")])

Menghitung proporsi dan jumlah setiap kategori dalam variabel target, untuk memastikan distribusi data seimbang. Jika adanya permasalahan data yang tidak seimbang, maka perlu dilakukan handling imbalance data.

round(prop.table(table(dataku$ASI_Eks)),4)
## 
##      0      1 
## 0.2884 0.7116

Handling Imbalancing Data

Jika dilihat dari proporsi kedua kategori di atas, tampak tidak cukup seimbang sehingga perlu membutuhkan pre-processing tambahan untuk menyeimbangkan proporsi antar dua kelas target variabel.

Variabel ASI Eksklusif merupakan variabel target dan yang bernilai 0 sangat sedikit (28,84%) dibandingkan dengan kelas mayoritas (71,16%). Oleh karena itu, selanjutnya dilakukan proses balancing dengan menggunakan teknik oversampling kelas minoritas, yaitu dengan menambahkan jumlah data pada kelas minoritas dengan cara menduplikasi data.

library(dplyr)

table(dataku$ASI_Eks)
## 
##    0    1 
## 1601 3950
# Oversampling
data_oversampled <- dataku %>%
  filter(ASI_Eks == 0) %>%
  slice(rep(1:n(), times = ceiling(sum(dataku$ASI_Eks == 1) / sum(dataku$ASI_Eks == 0)))) %>%
  bind_rows(filter(dataku, ASI_Eks == 1))

# Distribusi setelah oversampling
table(data_oversampled$ASI_Eks)
## 
##    0    1 
## 4803 3950

2 Model Regresi

probit <- glm(ASI_Eks~MISKIN+HH_FAMILY+R1401,data=data_oversampled, family = "binomial"(link="probit"))
summary(probit)
## 
## Call:
## glm(formula = ASI_Eks ~ MISKIN + HH_FAMILY + R1401, family = binomial(link = "probit"), 
##     data = data_oversampled)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.173715   0.068886   2.522 0.011677 *  
## MISKIN2     -0.146828   0.038546  -3.809 0.000139 ***
## HH_FAMILY2   0.247243   0.055507   4.454 8.42e-06 ***
## HH_FAMILY3   0.255738   0.057484   4.449 8.63e-06 ***
## R1401       -0.143574   0.008285 -17.330  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 12051  on 8752  degrees of freedom
## Residual deviance: 11703  on 8748  degrees of freedom
## AIC: 11713
## 
## Number of Fisher Scoring iterations: 4

Berdasarkan output di atas diperoleh bahwa status kemiskinan dan jumlah ART dalam rumah tangga berpengaruh secara nyata terhadap pemberian ASI Eksklusif bayi usia 0-5 bulan. Hal ini ditunjukkan dari nilai p-value yang diperoleh <0.05.

Pearson chi-square Test

Uji Pearson chi-square untuk goodness of fit bertujuan untuk menguji apakah data yang diamati cocok dengan distribusi yang diharapkan (model yang dibangun) (Hosmer & Lemeshow, 2000). Uji ini membandingkan frekuensi yang diamati dengan frekuensi yang diprediksi oleh model untuk melihat apakah ada perbedaan yang signifikan antara keduanya. 2

goodness_of_fit <- chisq.test(fitted(probit), data_oversampled$ASI_Eks)
print(goodness_of_fit)
## 
##  Pearson's Chi-squared test
## 
## data:  fitted(probit) and data_oversampled$ASI_Eks
## X-squared = 398.81, df = 35, p-value < 2.2e-16

Hasil uji chi-square menunjukkan nilai X-squared = 398.81 dengan p-value yang lebih kecil dari 0.05. Ini menunjukkan bahwa model probit yang dibangun cukup baik dalam mencocokkan data.

Interpretasi Hasil

probit$coefficient
## (Intercept)     MISKIN2  HH_FAMILY2  HH_FAMILY3       R1401 
##   0.1737155  -0.1468282   0.2472426   0.2557375  -0.1435741

Berikut adalah output koefisien model regresi probit:

Variabel Koefisien
Intercept 0.1737
MISKIN2 -0.1468
HH_FAMILY2 0.2472
HH_FAMILY3 0.2557
R1401 -0.1436

Koefisien regresi probit tidak secara langsung menunjukkan perubahan probabilitas tetapi menunjukkan efek variabel independen pada fungsi distribusi kumulatif normal standar.

1. Intercept (0.1737)

  • Nilai intercept mewakili fungsi tautan probit (z-score) saat semua variabel independen bernilai nol. Dalam konteks ini, nilai 0.1737 menunjukkan bahwa bayi yang lahir di rumah tangga dengan semua variabel independen pada nilai referensi memiliki probabilitas yang lebih besar untuk mendapatkan ASI eksklusif.

2. MISKIN2 (-0.1468)

  • Variabel MISKIN2 adalah indikator apakah keluarga tidak miskin (2 = tidak miskin, 1 = miskin).

  • Koefisien -0.1468 menunjukkan bahwa bayi dari rumah tangga tidak miskin memiliki probabilitas lebih kecil untuk memberikan ASI eksklusif dibandingkan bayi dari rumah tangga miskin.

  • Efek ini mengindikasikan bahwa keluarga miskin mungkin memiliki tingkat kepatuhan lebih tinggi terhadap pemberian ASI eksklusif, mungkin karena keterbatasan akses ke susu formula atau alasan lainnya.

3. HH_FAMILY2 (0.2472)

  • Variabel HH_FAMILY2 menunjukkan keluarga dengan 3–5 anggota rumah tangga (ART).

  • Koefisien 0.2472 menunjukkan bahwa bayi yang lahir dalam rumah tangga dengan 3–5 ART memiliki probabilitas lebih besar untuk memberikan ASI eksklusif dibandingkan dengan keluarga referensi (keluarga dengan kurang dari 3 ART).

  • Hal ini mungkin terjadi karena keluarga dengan ukuran sedang memiliki cukup sumber daya dan dukungan untuk pemberian ASI eksklusif.

4. HH_FAMILY3 (0.2557)

  • Variabel HH_FAMILY3 menunjukkan rumah tangga dengan 6 atau lebih anggota rumah tangga (ART).

  • Koefisien 0.2557 menunjukkan bahwa bayi dari rumah tangga dengan 6 atau lebih ART memiliki probabilitas lebih besar untuk mendapatkan ASI eksklusif dibandingkan dengan keluarga referensi.

5. R1401 (-0.1436)

  • Variabel R1401 adalah usia anak dalam bulan (0–59 bulan).

  • Koefisien -0.1436 menunjukkan bahwa setiap peningkatan 1 bulan usia anak mengurangi probabilitas mendapatkan ASI eksklusif.

  • Hal ini masuk akal karena pemberian ASI eksklusif biasanya hanya diberikan pada bayi usia 0–6 bulan, sehingga probabilitas menurun seiring bertambahnya usia anak.

3 Splitting Train-Test

Membagi data menjadi 80% untuk pelatihan (train) dan 20% untuk pengujian (test), dengan seed untuk memastikan hasil yang dapat direproduksi.

set.seed(785)
sample <- sample(nrow(data_oversampled),floor(nrow(data_oversampled)*0.8))
train <- data_oversampled[sample,]
dim(train)
## [1] 7002   17
test <- data_oversampled[-sample,]
dim(test)
## [1] 1751   17

4 Evaluasi Model

test$pred <- predict(probit, test, type="response")
test$ASI_Eks_pred <- ifelse(test$pred > 0.6, "1", "0")
test$ASI_Eks_pred <- as.factor(test$ASI_Eks_pred)
test$ASI_Eks_pred <- factor(test$ASI_Eks_pred, levels = levels(test$ASI_Eks))

Confusion Matrix

conf.mat <- caret::confusionMatrix(test$ASI_Eks_pred, test$ASI_Eks, positive = "0")
print(conf.mat)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 871 653
##          1  89 138
##                                           
##                Accuracy : 0.5762          
##                  95% CI : (0.5527, 0.5995)
##     No Information Rate : 0.5483          
##     P-Value [Acc > NIR] : 0.009827        
##                                           
##                   Kappa : 0.0872          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9073          
##             Specificity : 0.1745          
##          Pos Pred Value : 0.5715          
##          Neg Pred Value : 0.6079          
##              Prevalence : 0.5483          
##          Detection Rate : 0.4974          
##    Detection Prevalence : 0.8704          
##       Balanced Accuracy : 0.5409          
##                                           
##        'Positive' Class : 0               
## 
broom::tidy(conf.mat)
## # A tibble: 14 × 6
##    term                 class estimate conf.low conf.high   p.value
##    <chr>                <chr>    <dbl>    <dbl>     <dbl>     <dbl>
##  1 accuracy             <NA>    0.576     0.553     0.600  9.83e- 3
##  2 kappa                <NA>    0.0872   NA        NA     NA       
##  3 mcnemar              <NA>   NA        NA        NA      6.67e-95
##  4 sensitivity          0       0.907    NA        NA     NA       
##  5 specificity          0       0.174    NA        NA     NA       
##  6 pos_pred_value       0       0.572    NA        NA     NA       
##  7 neg_pred_value       0       0.608    NA        NA     NA       
##  8 precision            0       0.572    NA        NA     NA       
##  9 recall               0       0.907    NA        NA     NA       
## 10 f1                   0       0.701    NA        NA     NA       
## 11 prevalence           0       0.548    NA        NA     NA       
## 12 detection_rate       0       0.497    NA        NA     NA       
## 13 detection_prevalence 0       0.870    NA        NA     NA       
## 14 balanced_accuracy    0       0.541    NA        NA     NA

Hasil confusion matrix menunjukkan bahwa model regresi Probit yang diperoleh memiliki persentase pengamatan (aktual) yang diprediksi dengan benar oleh model sebesar 57.60 persen. Persentase dari kelas 0 (aktual) yang diprediksi dengan benar oleh model adalah sebesar 90.70 persen, sementara persentase dari kelas 1 (aktual) yang diprediksi dengan benar oleh model adalah sebesar 17.40 persen, dengan accuracy keseluruhan sebesar 57.60 persen.

Plot ROC Curve

ROC Curve (Receiver Operating Characteristic Curve) adalah grafik yang digunakan untuk mengevaluasi performa model klasifikasi, khususnya pada masalah dengan variabel target biner (0/1, Ya/Tidak, Positif/Negatif). Plot ini menunjukkan kemampuan model untuk membedakan antara kelas positif dan negatif pada berbagai ambang batas (threshold).

  • Kurva mendekati sudut kiri atas: Model sangat baik.

  • Kurva mendekati diagonal: Model mendekati tebakan acak.

library(pROC)
roc_obj <- roc(response = test$ASI_Eks, predictor = test$pred, levels = c("0", "1"))
plot(roc_obj, main = "ROC Curve", col = "blue", lwd = 2)

Menampilkan AUC

AUC (Area Under the Curve) adalah ukuran performa model klasifikasi berdasarkan kurva ROC.

Kurva ROC memplot hubungan antara:

  • True Positive Rate (TPR) atau sensitivitas.

  • False Positive Rate (FPR) (1 - spesifisitas).

Nilai AUC berkisar antara 0 dan 1:

  • 1.0: Model sempurna dalam membedakan kelas.

  • 0.5: Model tidak lebih baik daripada tebakan acak.

  • < 0.5: Model berkinerja lebih buruk daripada tebakan acak.

auc_value <- auc(roc_obj)
print(paste("AUC:", round(auc_value,4)))
## [1] "AUC: 0.5839"

Nilai AUC diperoleh sebesar 0.5839 yang menunjukkan model memiliki kemampuan moderat dalam membedakan antara kelas positif (1) dan negatif (0).

Dalam konteks ini, jika memilih dua pengamatan secara acak, satu dari kelas positif (1) dan satu dari kelas negatif (0), model memiliki probabilitas sekitar 58.39% untuk mengklasifikasikan pengamatan positif dengan probabilitas lebih tinggi daripada pengamatan negatif.

5 References

[1] Agresti A. (2002). Categorical Data Analysis. New York : John Willey and Sons.

[2] Hosmer, D. W., & Lemeshow, S. (2000). Applied Logistic Regression. Wiley-Interscience.

Direktorat Statistik Kesejahteraan Rakyat, BPS,

  1. Yang bisa ditampilkan dalam plot bar adalah variabel yang bertipe factor.↩︎

  2. Dalam konteks regresi, ini digunakan untuk menilai sejauh mana model yang dibangun dapat menggambarkan data yang sesungguhnya.↩︎