Regresi Logistik Multinomial Preferensi KB

Deskripsi Data

Data Preferensi KB

Daftar Variabel

No	Var	Keterangan
X1	Umur	Numeric
X2	Pendidikan	0-SMP kebawah 1-SMA 2-PT
X3	Status Pekerjaan	0-Tidak bekerja 1-Bekerja
X4	Status Melahirkan	0-Tidak pernah 1-Pernah
Y	Preferensi KB	1-Pil 2-Suntik 3-IUD

Package

library(haven)
library(nnet)
library(readxl)
library(caret)

## Loading required package: ggplot2

## Loading required package: lattice

library(arsenal)

1. Input Data

data <- read_xlsx("DATA MULTINOM.xlsx", sheet = "Preferensi KB")
data

## # A tibble: 9,774 × 6
##    Status.KRT   Umur Pendidikan    Status.Kerja  Status.Melahirkan Preferensi.KB
##    <chr>       <dbl> <chr>         <chr>         <chr>             <chr>        
##  1 0-Bukan KRT    24 0-SMP kebawah 1-Bekerja     1-Pernah          2-Suntik     
##  2 0-Bukan KRT    45 0-SMP kebawah 1-Bekerja     1-Pernah          3-IUD        
##  3 0-Bukan KRT    48 0-SMP kebawah 0-Tidak Beke… 1-Pernah          3-IUD        
##  4 0-Bukan KRT    51 0-SMP kebawah 1-Bekerja     1-Pernah          3-IUD        
##  5 0-Bukan KRT    38 1-SMA         1-Bekerja     1-Pernah          3-IUD        
##  6 0-Bukan KRT    24 1-SMA         0-Tidak Beke… 1-Pernah          2-Suntik     
##  7 0-Bukan KRT    47 0-SMP kebawah 1-Bekerja     1-Pernah          2-Suntik     
##  8 0-Bukan KRT    22 0-SMP kebawah 0-Tidak Beke… 1-Pernah          1-Pil        
##  9 0-Bukan KRT    35 1-SMA         0-Tidak Beke… 1-Pernah          2-Suntik     
## 10 0-Bukan KRT    42 0-SMP kebawah 1-Bekerja     1-Pernah          2-Suntik     
## # ℹ 9,764 more rows

1.1. Struktur data

str(data)

## tibble [9,774 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Status.KRT       : chr [1:9774] "0-Bukan KRT" "0-Bukan KRT" "0-Bukan KRT" "0-Bukan KRT" ...
##  $ Umur             : num [1:9774] 24 45 48 51 38 24 47 22 35 42 ...
##  $ Pendidikan       : chr [1:9774] "0-SMP kebawah" "0-SMP kebawah" "0-SMP kebawah" "0-SMP kebawah" ...
##  $ Status.Kerja     : chr [1:9774] "1-Bekerja" "1-Bekerja" "0-Tidak Bekerja" "1-Bekerja" ...
##  $ Status.Melahirkan: chr [1:9774] "1-Pernah" "1-Pernah" "1-Pernah" "1-Pernah" ...
##  $ Preferensi.KB    : chr [1:9774] "2-Suntik" "3-IUD" "3-IUD" "3-IUD" ...

1.2. Merubah ke factor

data$Status.KRT <- as.factor(data$Status.KRT)

data$Pendidikan <- as.factor(data$Pendidikan)
data$Status.Kerja <- as.factor(data$Status.Kerja)
data$Status.Melahirkan <- as.factor(data$Status.Melahirkan)
data$Preferensi.KB <- as.factor(data$Preferensi.KB)


str(data)

## tibble [9,774 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Status.KRT       : Factor w/ 2 levels "0-Bukan KRT",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Umur             : num [1:9774] 24 45 48 51 38 24 47 22 35 42 ...
##  $ Pendidikan       : Factor w/ 3 levels "0-SMP kebawah",..: 1 1 1 1 2 2 1 1 2 1 ...
##  $ Status.Kerja     : Factor w/ 2 levels "0-Tidak Bekerja",..: 2 2 1 2 2 1 2 1 1 2 ...
##  $ Status.Melahirkan: Factor w/ 2 levels "0-Tidak pernah",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Preferensi.KB    : Factor w/ 3 levels "1-Pil","2-Suntik",..: 2 3 3 3 3 2 2 1 2 2 ...

2. EDA

2.1. Plot

boxplot(Umur~ Preferensi.KB, data = data,
        main = "Box Plot Umur dan Preferensi KB",
        col = c("lightcyan","lavender", "mistyrose")
        )

# Preferensi KB
barplot(table(data$Preferensi.KB),
        main = "Bar Plot Preferensi KB",
        xlab = "Preferensi KB",
        ylab = "Jumlah",
        col = c("lightcyan","lavender", "mistyrose")
        )

# Preferensi KB dan Pendidikan
barplot(table(data$Preferensi.KB, data$Pendidikan),
        main = "Bar Plot Pendidikan dan Preferensi KB",
        xlab = "Pendidikan",
        ylab = "Jumlah",
        col = c("lightcyan","lavender", "mistyrose"),
        beside = TRUE
        )
legend("topright",
       c("Pil","Suntik", "IUD"),
       fill = c("lightcyan","lavender", "mistyrose")
       )

# Preferensi KB dan status kerja
barplot(table(data$Preferensi.KB, data$Status.Kerja),
        main = "Bar Plot Status Pekerjaan dan Preferensi KB",
        xlab = "Pekerjaan",
        ylab = "Jumlah",
        col = c("lightcyan","lavender", "mistyrose"),
        beside = TRUE
        )
legend("topright",
       c("Pil","Suntik", "IUD"),
       fill = c("lightcyan","lavender", "mistyrose")
       )

2.2. Summary Table

tab<-tableby(Preferensi.KB~ .,data = data)
summary(tab,text=TRUE)

## 
## 
## |                   | 1-Pil (N=2441)  | 2-Suntik (N=5264) | 3-IUD (N=2069)  | Total (N=9774)  | p value|
## |:------------------|:---------------:|:-----------------:|:---------------:|:---------------:|-------:|
## |Status.KRT         |                 |                   |                 |                 |   0.009|
## |-  0-Bukan KRT     |  2428 (99.5%)   |   5233 (99.4%)    |  2044 (98.8%)   |  9705 (99.3%)   |        |
## |-  1-KRT           |    13 (0.5%)    |     31 (0.6%)     |    25 (1.2%)    |    69 (0.7%)    |        |
## |Umur               |                 |                   |                 |                 | < 0.001|
## |-  Mean (SD)       | 37.034 (8.137)  |  33.978 (8.232)   | 38.076 (7.612)  | 35.609 (8.278)  |        |
## |-  Range           | 15.000 - 54.000 |  14.000 - 54.000  | 16.000 - 54.000 | 14.000 - 54.000 |        |
## |Pendidikan         |                 |                   |                 |                 | < 0.001|
## |-  0-SMP kebawah   |  1745 (71.5%)   |   3701 (70.3%)    |   997 (48.2%)   |  6443 (65.9%)   |        |
## |-  1-SMA           |   584 (23.9%)   |   1313 (24.9%)    |   671 (32.4%)   |  2568 (26.3%)   |        |
## |-  2-PT            |   112 (4.6%)    |    250 (4.7%)     |   401 (19.4%)   |   763 (7.8%)    |        |
## |Status.Kerja       |                 |                   |                 |                 | < 0.001|
## |-  0-Tidak Bekerja |  1424 (58.3%)   |   3483 (66.2%)    |  1152 (55.7%)   |  6059 (62.0%)   |        |
## |-  1-Bekerja       |  1017 (41.7%)   |   1781 (33.8%)    |   917 (44.3%)   |  3715 (38.0%)   |        |
## |Status.Melahirkan  |                 |                   |                 |                 | < 0.001|
## |-  0-Tidak pernah  |    44 (1.8%)    |     41 (0.8%)     |    5 (0.2%)     |    90 (0.9%)    |        |
## |-  1-Pernah        |  2397 (98.2%)   |   5223 (99.2%)    |  2064 (99.8%)   |  9684 (99.1%)   |        |

3. Partition data

set.seed(1234)

acak <- createDataPartition(data$Preferensi.KB, p=0.8, list=FALSE)
data_train <- data[acak,]
data_test <- data[-acak,]

summary(data_train)

##        Status.KRT        Umur               Pendidikan            Status.Kerja 
##  0-Bukan KRT:7762   Min.   :15.00   0-SMP kebawah:5143   0-Tidak Bekerja:4835  
##  1-KRT      :  59   1st Qu.:29.00   1-SMA        :2075   1-Bekerja      :2986  
##                     Median :36.00   2-PT         : 603                         
##                     Mean   :35.63                                              
##                     3rd Qu.:42.00                                              
##                     Max.   :54.00                                              
##       Status.Melahirkan  Preferensi.KB 
##  0-Tidak pernah:  74    1-Pil   :1953  
##  1-Pernah      :7747    2-Suntik:4212  
##                         3-IUD   :1656  
##                                        
##                                        
##

summary(data_test)

##        Status.KRT        Umur               Pendidikan            Status.Kerja 
##  0-Bukan KRT:1943   Min.   :14.00   0-SMP kebawah:1300   0-Tidak Bekerja:1224  
##  1-KRT      :  10   1st Qu.:29.00   1-SMA        : 493   1-Bekerja      : 729  
##                     Median :36.00   2-PT         : 160                         
##                     Mean   :35.52                                              
##                     3rd Qu.:42.00                                              
##                     Max.   :54.00                                              
##       Status.Melahirkan  Preferensi.KB 
##  0-Tidak pernah:  16    1-Pil   : 488  
##  1-Pernah      :1937    2-Suntik:1052  
##                         3-IUD   : 413  
##                                        
##                                        
##

4. Create Model

model <- multinom(Preferensi.KB~ Umur+Pendidikan+Status.Kerja, data = data_train)

## # weights:  18 (10 variable)
## initial  value 8592.246710 
## iter  10 value 7599.353956
## final  value 7476.893005 
## converged

summary(model)

## Call:
## multinom(formula = Preferensi.KB ~ Umur + Pendidikan + Status.Kerja, 
##     data = data_train)
## 
## Coefficients:
##          (Intercept)        Umur Pendidikan1-SMA Pendidikan2-PT
## 2-Suntik    2.410224 -0.04368472      0.03056896      0.0733296
## 3-IUD      -1.298157  0.02061109      0.70793298      1.8381077
##          Status.Kerja1-Bekerja
## 2-Suntik           -0.26418542
## 3-IUD              -0.08812783
## 
## Std. Errors:
##          (Intercept)        Umur Pendidikan1-SMA Pendidikan2-PT
## 2-Suntik   0.1295495 0.003463652      0.06516110      0.1329431
## 3-IUD      0.1713646 0.004357311      0.07759096      0.1286837
##          Status.Kerja1-Bekerja
## 2-Suntik            0.05769147
## 3-IUD               0.07060176
## 
## Residual Deviance: 14953.79 
## AIC: 14973.79

4.1. Uji multicol

library(car)

vif(model)

## Warning in vif.default(model): No intercept: vifs may not be sensible.

##                   GVIF Df GVIF^(1/(2*Df))
## Umur         29.562900  1        5.437178
## Pendidikan    4.594105  2        1.464031
## Status.Kerja  2.337103  1        1.528759

4.2. Uji simultan

library(lmtest)
lrtest(model)

## # weights:  6 (2 variable)
## initial  value 8592.246710 
## final  value 7887.168974 
## converged

## Likelihood ratio test
## 
## Model 1: Preferensi.KB ~ Umur + Pendidikan + Status.Kerja
## Model 2: Preferensi.KB ~ 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1  10 -7476.9                         
## 2   2 -7887.2 -8 820.55  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.3. Gof Test

library(pscl)
pR2(model)

## fitting null model for pseudo-r2
## # weights:  6 (2 variable)
## initial  value 8592.246710 
## final  value 7887.168974 
## converged

##           llh       llhNull            G2      McFadden          r2ML 
## -7.476893e+03 -7.887169e+03  8.205519e+02  5.201815e-02  9.960030e-02 
##          r2CU 
##  1.148878e-01

4.4. Hitung p-value menggunakan z test

z <- summary(model)$coefficients/summary(model)$standard.errors
p <- (1-pnorm(abs(z), 0, 1))*2
data.frame(p)

##          X.Intercept.         Umur Pendidikan1.SMA Pendidikan2.PT
## 2-Suntik 0.000000e+00 0.000000e+00       0.6389774      0.5812317
## 3-IUD    3.574918e-14 2.242642e-06       0.0000000      0.0000000
##          Status.Kerja1.Bekerja
## 2-Suntik          4.665777e-06
## 3-IUD             2.119438e-01

5. Prediksi pada data test

predict_prob = predict(model, data_test, type = "prob")
head(predict_prob,10)

##        1-Pil  2-Suntik     3-IUD
## 1  0.3223394 0.4409638 0.2366968
## 2  0.2969277 0.5503545 0.1527178
## 3  0.2175926 0.5897773 0.1926301
## 4  0.2640707 0.4426615 0.2932678
## 5  0.2488155 0.4967263 0.2544582
## 6  0.2691809 0.5699729 0.1608461
## 7  0.3077062 0.4798903 0.2124034
## 8  0.2396153 0.6312257 0.1291589
## 9  0.1614858 0.2478224 0.5906918
## 10 0.1425234 0.2391876 0.6182890

head(data.frame(predict_prob, data_test$Preferensi.KB),10)

##       X1.Pil X2.Suntik    X3.IUD data_test.Preferensi.KB
## 1  0.3223394 0.4409638 0.2366968                   3-IUD
## 2  0.2969277 0.5503545 0.1527178                2-Suntik
## 3  0.2175926 0.5897773 0.1926301                2-Suntik
## 4  0.2640707 0.4426615 0.2932678                   3-IUD
## 5  0.2488155 0.4967263 0.2544582                2-Suntik
## 6  0.2691809 0.5699729 0.1608461                2-Suntik
## 7  0.3077062 0.4798903 0.2124034                   1-Pil
## 8  0.2396153 0.6312257 0.1291589                2-Suntik
## 9  0.1614858 0.2478224 0.5906918                2-Suntik
## 10 0.1425234 0.2391876 0.6182890                   3-IUD

5.1. Confusion Matrix

prediksi.test <- predict(model, data_test,type = "class")
data_test$Preferensi.KB<-as.factor(data_test$Preferensi.KB)
confusionMatrix(as.factor(prediksi.test),
                data_test$Preferensi.KB)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 1-Pil 2-Suntik 3-IUD
##   1-Pil        8        9     5
##   2-Suntik   450      990   313
##   3-IUD       30       53    95
## 
## Overall Statistics
##                                           
##                Accuracy : 0.5597          
##                  95% CI : (0.5373, 0.5818)
##     No Information Rate : 0.5387          
##     P-Value [Acc > NIR] : 0.03289         
##                                           
##                   Kappa : 0.1094          
##                                           
##  Mcnemar's Test P-Value : < 2e-16         
## 
## Statistics by Class:
## 
##                      Class: 1-Pil Class: 2-Suntik Class: 3-IUD
## Sensitivity              0.016393          0.9411      0.23002
## Specificity              0.990444          0.1532      0.94610
## Pos Pred Value           0.363636          0.5647      0.53371
## Neg Pred Value           0.751424          0.6900      0.82085
## Prevalence               0.249872          0.5387      0.21147
## Detection Rate           0.004096          0.5069      0.04864
## Detection Prevalence     0.011265          0.8976      0.09114
## Balanced Accuracy        0.503419          0.5471      0.58806

6. Interpretasi koefisien regresi

data.frame(summary(model)$coefficients)

##          X.Intercept.        Umur Pendidikan1.SMA Pendidikan2.PT
## 2-Suntik     2.410224 -0.04368472      0.03056896      0.0733296
## 3-IUD       -1.298157  0.02061109      0.70793298      1.8381077
##          Status.Kerja1.Bekerja
## 2-Suntik           -0.26418542
## 3-IUD              -0.08812783

# odds = exponensial parameter model
data.frame(exp(summary(model)$coefficients))

##          X.Intercept.      Umur Pendidikan1.SMA Pendidikan2.PT
## 2-Suntik   11.1364601 0.9572557        1.031041       1.076085
## 3-IUD       0.2730345 1.0208250        2.029791       6.284634
##          Status.Kerja1.Bekerja
## 2-Suntik             0.7678312
## 3-IUD                0.9156438

KB Suntik terhadap Pil

Variabel	Odds	Keterangan
Intercept	11.1364601	Responden dengan pendidikan SMP dan tidak bekerja, odds memilih KB suntik dibanding pil adalah 11.136.
Umur	0.9572557	Setiap penambahan umur 1 tahun, odds responden memilih KB suntik adalah 0.957 kali dibanding memilih KB pil.
Pendidikan (SMA)	1.031041	Responden dengan pendidikan SMA, odds responden memilih KB suntik adalah 1.031 kalinya dibanding odds responden memilih KB pil.
Pendidikan (PT)	1.076085	Responden dengan pendidikan PT, odds responden memilih KB suntik adalah 1.076 kalinya dibanding odds responden memilih KB pil.
Status pekerjaan (Bekerja)	0.7678312	Responden dengan status bekerja, odds responden memilih KB suntik adalah 0.7678 kalinya dibanding odds responden memilih KB pil.

KB IUD terhadap Pil

Variabel	Odds	Keterangan
Intercept	0.2730345	Responden dengan pendidikan SMP dan tidak bekerja, odds memilih IUD suntik dibanding pil adalah 0.2730345.
Umur	1.0208250	Setiap penambahan umur 1 tahun, odds responden memilih KB IUDadalah 1.0208 kali dibanding memilih KB pil.
Pendidikan (SMA)	2.029791	Responden dengan pendidikan SMA, odds responden memilih KB IUD adalah 2.0297 kalinya dibanding odds responden memilih KB pil.
Pendidikan (PT)	6.284634	Responden dengan pendidikan PT, odds responden memilih KB IUD adalah 6.2846 kalinya dibanding odds responden memilih KB pil.
Status pekerjaan (Bekerja)	0.9156438	Responden dengan status bekerja, odds responden memilih KB IUD adalah 0.9156 kalinya dibanding odds responden memilih KB pil.