Regresi Logistik Ordinal Pendapatan

Package

library(ggplot2)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(gmodels)
library(ggmosaic)
library(corrplot)

## corrplot 0.92 loaded

library(caret)

## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift

library(rpart)
library(rpart.plot)
library(C50)
library(gmodels)
library(tidyr)
library(dplyr)
library(knitr)
library(rsample)
library(readxl)
library(MASS)

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select

1. Input Data

1.1. Input

data_pendapatan<-read_xlsx("Data Latihan.xlsx", sheet = "Pendapatan")
head(data_pendapatan)

## # A tibble: 6 × 7
##   Pendidikan Klasifikasi.Daerah Jenis.Kelamin Status.KRT  Usia Status.Perkawinan
##        <dbl>              <dbl>         <dbl>      <dbl> <dbl>             <dbl>
## 1          0                  0             1          1    49                 1
## 2          0                  0             1          1    47                 1
## 3          0                  0             1          1    61                 1
## 4          1                  0             0          0    23                 0
## 5          0                  0             1          1    27                 1
## 6          0                  0             1          1    50                 1
## # ℹ 1 more variable: Pendapatan <dbl>

1.2. Keterangan

Variabel	Keterangan	Nilai
Y	Tingkat Pendapatan	Rendah Menengah kebawah Menengah keatas Tinggi
X1	Tingkat Pendidikan	SMA kebawah PT
X2	Klasifikasi desa kota	Kota Desa
X3	Jenis kelamin	Perempuan Laki-laki
X4	Status KRT	Bukan KRT KRT
X5	Usia	Numeric
X6	Status Perkawinan	Belum kawin Kawin

Sumber data dari SAKERNAS 2020

1.3. Merubah variabel ke factor

str(data_pendapatan)

## tibble [1,000 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Pendidikan        : num [1:1000] 0 0 0 1 0 0 0 0 0 0 ...
##  $ Klasifikasi.Daerah: num [1:1000] 0 0 0 0 0 0 1 1 1 1 ...
##  $ Jenis.Kelamin     : num [1:1000] 1 1 1 0 1 1 1 1 0 1 ...
##  $ Status.KRT        : num [1:1000] 1 1 1 0 1 1 1 1 0 0 ...
##  $ Usia              : num [1:1000] 49 47 61 23 27 50 46 49 45 43 ...
##  $ Status.Perkawinan : num [1:1000] 1 1 1 0 1 1 1 1 1 0 ...
##  $ Pendapatan        : num [1:1000] 2 2 1 2 2 2 4 1 1 1 ...

data_pendapatan$Pendapatan <- as.factor(data_pendapatan$Pendapatan)
data_pendapatan$Pendidikan <- as.factor(data_pendapatan$Pendidikan)
data_pendapatan$Klasifikasi.Daerah <- as.factor(data_pendapatan$Klasifikasi.Daerah)
data_pendapatan$Jenis.Kelamin <- as.factor(data_pendapatan$Jenis.Kelamin)
data_pendapatan$Status.KRT <- as.factor(data_pendapatan$Status.KRT)
data_pendapatan$Status.Perkawinan <- as.factor(data_pendapatan$Status.Perkawinan)

data_pendapatan

## # A tibble: 1,000 × 7
##    Pendidikan Klasifikasi.Daerah Jenis.Kelamin Status.KRT  Usia
##    <fct>      <fct>              <fct>         <fct>      <dbl>
##  1 0          0                  1             1             49
##  2 0          0                  1             1             47
##  3 0          0                  1             1             61
##  4 1          0                  0             0             23
##  5 0          0                  1             1             27
##  6 0          0                  1             1             50
##  7 0          1                  1             1             46
##  8 0          1                  1             1             49
##  9 0          1                  0             0             45
## 10 0          1                  1             0             43
## # ℹ 990 more rows
## # ℹ 2 more variables: Status.Perkawinan <fct>, Pendapatan <fct>

str(data_pendapatan)

## tibble [1,000 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Pendidikan        : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 1 1 1 ...
##  $ Klasifikasi.Daerah: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 2 2 2 ...
##  $ Jenis.Kelamin     : Factor w/ 2 levels "0","1": 2 2 2 1 2 2 2 2 1 2 ...
##  $ Status.KRT        : Factor w/ 2 levels "0","1": 2 2 2 1 2 2 2 2 1 1 ...
##  $ Usia              : num [1:1000] 49 47 61 23 27 50 46 49 45 43 ...
##  $ Status.Perkawinan : Factor w/ 2 levels "0","1": 2 2 2 1 2 2 2 2 2 1 ...
##  $ Pendapatan        : Factor w/ 4 levels "1","2","3","4": 2 2 1 2 2 2 4 1 1 1 ...

summary(data_pendapatan)

##  Pendidikan Klasifikasi.Daerah Jenis.Kelamin Status.KRT      Usia      
##  0:874      0:196              0:309         0:407      Min.   :14.00  
##  1:126      1:804              1:691         1:593      1st Qu.:30.00  
##                                                         Median :39.00  
##                                                         Mean   :39.89  
##                                                         3rd Qu.:49.00  
##                                                         Max.   :80.00  
##  Status.Perkawinan Pendapatan
##  0:262             1:279     
##  1:738             2:568     
##                    3:133     
##                    4: 20     
##                              
##

library(arsenal)

## 
## Attaching package: 'arsenal'

## The following object is masked from 'package:lubridate':
## 
##     is.Date

tab<-tableby(Pendapatan~ .,data = data_pendapatan)
summary(tab,text=TRUE)

## 
## 
## |                   |    1 (N=279)    |    2 (N=568)    |    3 (N=133)    |    4 (N=20)     | Total (N=1000)  | p value|
## |:------------------|:---------------:|:---------------:|:---------------:|:---------------:|:---------------:|-------:|
## |Pendidikan         |                 |                 |                 |                 |                 | < 0.001|
## |-  0               |   268 (96.1%)   |   524 (92.3%)   |   74 (55.6%)    |    8 (40.0%)    |   874 (87.4%)   |        |
## |-  1               |    11 (3.9%)    |    44 (7.7%)    |   59 (44.4%)    |   12 (60.0%)    |   126 (12.6%)   |        |
## |Klasifikasi.Daerah |                 |                 |                 |                 |                 | < 0.001|
## |-  0               |   66 (23.7%)    |   128 (22.5%)   |    2 (1.5%)     |    0 (0.0%)     |   196 (19.6%)   |        |
## |-  1               |   213 (76.3%)   |   440 (77.5%)   |   131 (98.5%)   |   20 (100.0%)   |   804 (80.4%)   |        |
## |Jenis.Kelamin      |                 |                 |                 |                 |                 | < 0.001|
## |-  0               |   132 (47.3%)   |   128 (22.5%)   |   45 (33.8%)    |    4 (20.0%)    |   309 (30.9%)   |        |
## |-  1               |   147 (52.7%)   |   440 (77.5%)   |   88 (66.2%)    |   16 (80.0%)    |   691 (69.1%)   |        |
## |Status.KRT         |                 |                 |                 |                 |                 | < 0.001|
## |-  0               |   141 (50.5%)   |   208 (36.6%)   |   54 (40.6%)    |    4 (20.0%)    |   407 (40.7%)   |        |
## |-  1               |   138 (49.5%)   |   360 (63.4%)   |   79 (59.4%)    |   16 (80.0%)    |   593 (59.3%)   |        |
## |Usia               |                 |                 |                 |                 |                 | < 0.001|
## |-  Mean (SD)       | 42.957 (14.304) | 38.165 (12.304) | 39.910 (10.340) | 46.000 (5.554)  | 39.891 (12.753) |        |
## |-  Range           | 15.000 - 80.000 | 14.000 - 78.000 | 19.000 - 66.000 | 37.000 - 56.000 | 14.000 - 80.000 |        |
## |Status.Perkawinan  |                 |                 |                 |                 |                 |   0.007|
## |-  0               |   91 (32.6%)    |   138 (24.3%)   |   32 (24.1%)    |    1 (5.0%)     |   262 (26.2%)   |        |
## |-  1               |   188 (67.4%)   |   430 (75.7%)   |   101 (75.9%)   |   19 (95.0%)    |   738 (73.8%)   |        |

2. Pembentukan Model

2.1. Partisi data

set.seed(123)
acak <- createDataPartition(data_pendapatan$Pendapatan, p=0.8, list=FALSE)
data_train <- data_pendapatan[acak,]
data_test <- data_pendapatan[-acak,]

summary(data_train)

##  Pendidikan Klasifikasi.Daerah Jenis.Kelamin Status.KRT      Usia      
##  0:700      0:164              0:243         0:329      Min.   :15.00  
##  1:102      1:638              1:559         1:473      1st Qu.:30.00  
##                                                         Median :39.00  
##                                                         Mean   :39.57  
##                                                         3rd Qu.:48.00  
##                                                         Max.   :80.00  
##  Status.Perkawinan Pendapatan
##  0:219             1:224     
##  1:583             2:455     
##                    3:107     
##                    4: 16     
##                              
##

summary(data_test)

##  Pendidikan Klasifikasi.Daerah Jenis.Kelamin Status.KRT      Usia      
##  0:174      0: 32              0: 66         0: 78      Min.   :14.00  
##  1: 24      1:166              1:132         1:120      1st Qu.:32.00  
##                                                         Median :40.00  
##                                                         Mean   :41.18  
##                                                         3rd Qu.:49.75  
##                                                         Max.   :78.00  
##  Status.Perkawinan Pendapatan
##  0: 43             1: 55     
##  1:155             2:113     
##                    3: 26     
##                    4:  4     
##                              
##

2.2. Model

model_pendapatan<-polr(Pendapatan~., method = "logistic", data = data_train, Hess = T)
summary(model_pendapatan)

## Call:
## polr(formula = Pendapatan ~ ., data = data_train, Hess = T, method = "logistic")
## 
## Coefficients:
##                        Value Std. Error t value
## Pendidikan1          2.52138   0.237066  10.636
## Klasifikasi.Daerah1  0.56414   0.175636   3.212
## Jenis.Kelamin1       0.54822   0.199061   2.754
## Status.KRT1          0.64226   0.212620   3.021
## Usia                -0.03221   0.006865  -4.692
## Status.Perkawinan1   0.53006   0.180398   2.938
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -0.5217  0.3040    -1.7161
## 2|3  2.6805  0.3220     8.3249
## 3|4  5.1927  0.4173    12.4437
## 
## Residual Deviance: 1454.495 
## AIC: 1472.495

2.2.1. Uji asumsi multikolinieritas

#uji asumsi multikolinieritas

library(car)
vif(model_pendapatan)

##         Pendidikan Klasifikasi.Daerah      Jenis.Kelamin         Status.KRT 
##           1.045494           1.026166           1.637565           2.146962 
##               Usia  Status.Perkawinan 
##           1.467146           1.279745

2.2.2. Uji GoF

#Uji GoF

library(generalhoslem)
lipsitz.test(model_pendapatan)

## 
##  Lipsitz goodness of fit test for ordinal response models
## 
## data:  formula:  Pendapatan ~ Pendidikan + Klasifikasi.Daerah + Jenis.Kelamin + formula:      Status.KRT + Usia + Status.Perkawinan
## LR statistic = 25.252, df = 9, p-value = 0.002705

2.2.3. Simultant Test

library(lmtest) 
lrtest(model_pendapatan)

## Likelihood ratio test
## 
## Model 1: Pendapatan ~ Pendidikan + Klasifikasi.Daerah + Jenis.Kelamin + 
##     Status.KRT + Usia + Status.Perkawinan
## Model 2: Pendapatan ~ 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   9 -727.25                         
## 2   3 -821.76 -6 189.03  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.2.4. Partial Test

coefmodel<-c(-model_pendapatan$coefficients)
koefisien<-coef(summary(model_pendapatan))

#menghitung pvalue 
p <- pnorm(abs(koefisien[,"t value"]), lower.tail = FALSE)*2
(ctabel<-cbind(round(koefisien,4), "pvalue"=round(p,4)))

##                       Value Std. Error t value pvalue
## Pendidikan1          2.5214     0.2371 10.6358 0.0000
## Klasifikasi.Daerah1  0.5641     0.1756  3.2120 0.0013
## Jenis.Kelamin1       0.5482     0.1991  2.7540 0.0059
## Status.KRT1          0.6423     0.2126  3.0207 0.0025
## Usia                -0.0322     0.0069 -4.6924 0.0000
## Status.Perkawinan1   0.5301     0.1804  2.9383 0.0033
## 1|2                 -0.5217     0.3040 -1.7161 0.0862
## 2|3                  2.6805     0.3220  8.3249 0.0000
## 3|4                  5.1927     0.4173 12.4437 0.0000

3. Prediction of Data Test

#1
predict_prob = predict(model_pendapatan, data_test, type = "prob")
head(predict_prob,10)

##             1         2          3           4
## 1  0.33989193 0.5869007 0.06684284 0.006364482
## 2  0.09094504 0.6200168 0.25712374 0.031914397
## 3  0.34715626 0.5817922 0.06488758 0.006163973
## 4  0.22654704 0.6515250 0.11079339 0.011134596
## 5  0.45852474 0.4956463 0.04194932 0.003879622
## 6  0.27736764 0.6268207 0.08729227 0.008519372
## 7  0.45053826 0.5022034 0.04325224 0.004006122
## 8  0.42708974 0.5211741 0.04733153 0.004404657
## 9  0.16596916 0.6643290 0.15339845 0.016303380
## 10 0.47871059 0.4788780 0.03883282 0.003578580

head(data.frame(predict_prob, data_test$Pendapatan),20)

##            X1        X2         X3          X4 data_test.Pendapatan
## 1  0.33989193 0.5869007 0.06684284 0.006364482                    2
## 2  0.09094504 0.6200168 0.25712374 0.031914397                    2
## 3  0.34715626 0.5817922 0.06488758 0.006163973                    2
## 4  0.22654704 0.6515250 0.11079339 0.011134596                    1
## 5  0.45852474 0.4956463 0.04194932 0.003879622                    1
## 6  0.27736764 0.6268207 0.08729227 0.008519372                    3
## 7  0.45053826 0.5022034 0.04325224 0.004006122                    2
## 8  0.42708974 0.5211741 0.04733153 0.004404657                    2
## 9  0.16596916 0.6643290 0.15339845 0.016303380                    1
## 10 0.47871059 0.4788780 0.03883282 0.003578580                    1
## 11 0.41888599 0.5277037 0.04885576 0.004554542                    3
## 12 0.17047630 0.6643126 0.14941636 0.015794731                    2
## 13 0.19957107 0.6601811 0.12719275 0.013055090                    3
## 14 0.23224118 0.6492378 0.10773556 0.010785439                    1
## 15 0.49064332 0.4688438 0.03710068 0.003412192                    1
## 16 0.16155797 0.6641529 0.15746103 0.016828130                    2
## 17 0.17978138 0.6637024 0.14169240 0.014823835                    2
## 18 0.06370555 0.5621674 0.32789533 0.046231692                    3
## 19 0.45053826 0.5022034 0.04325224 0.004006122                    2
## 20 0.02126057 0.3268849 0.52003914 0.131815352                    2

3.1. Confussion Matrix of Model

#1
prediksi.test <- predict(model_pendapatan, data_test,type = "class")
data_test$Pendapatan<-as.factor(data_test$Pendapatan)
confusionMatrix(as.factor(prediksi.test), 
                data_test$Pendapatan, positive = "1")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   1   2   3   4
##          1  11   2   0   0
##          2  43 107  23   2
##          3   1   4   3   2
##          4   0   0   0   0
## 
## Overall Statistics
##                                           
##                Accuracy : 0.6111          
##                  95% CI : (0.5394, 0.6794)
##     No Information Rate : 0.5707          
##     P-Value [Acc > NIR] : 0.1406          
##                                           
##                   Kappa : 0.1738          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3 Class: 4
## Sensitivity           0.20000   0.9469  0.11538   0.0000
## Specificity           0.98601   0.2000  0.95930   1.0000
## Pos Pred Value        0.84615   0.6114  0.30000      NaN
## Neg Pred Value        0.76216   0.7391  0.87766   0.9798
## Prevalence            0.27778   0.5707  0.13131   0.0202
## Detection Rate        0.05556   0.5404  0.01515   0.0000
## Detection Prevalence  0.06566   0.8838  0.05051   0.0000
## Balanced Accuracy     0.59301   0.5735  0.53734   0.5000

Nilai akurasi sebesar 0.6111, artinya ketepatan model dalam melakukan prediksi sebesar 61.11%

3.2. Odds ratio

Untuk menginterpretasikan model, maka dapat dilakukan dengan menghitung nilai odds rasio terlebih dahulu.

#Odds Ratio
exp(coef(model_pendapatan))

##         Pendidikan1 Klasifikasi.Daerah1      Jenis.Kelamin1         Status.KRT1 
##          12.4456991           1.7579291           1.7301692           1.9007755 
##                Usia  Status.Perkawinan1 
##           0.9683003           1.6990306

Dari output tersebut, diperoleh informasi sebagai berikut.

Variabel	Keterangan	Odds rasio	Interpretasi
X1	Tingkat Pendidikan	12.4	Responden dengan pendidikan perguruan tinggi memiliki peluang 12.4 kali lebih besar untuk semakin mendapatkan pendapatan yang lebih tinggi dibanding lulusan SMA ke bawah.
X2	Klasifikasi desa kota	1.76	Responden yang tinggal di daerah desa justru memiliki peluang 1.76 kali lebih besar untuk menghasilkan pendapatan yang lebih besar dibanding penduduk yang tinggal di kota
X3	Jenis kelamin	1.73	Responden laki-laki memiliki peluang 1.73 kali lebih besar dibanding responden perempuan dalam mendapatkan pendapatan yang lebih besar.
X4	Status KRT	1.9	Responden dengan status kepala rumah tangga memiliki peluang 1.9 kali lebih besar dibanding responden dengan status bukan KRT dalam mendapatkan pendapatan yang lebih besar.
X5	Usia	0.97	Umur responden bertambah 1 tahun, maka kecenderungan ia untuk mendapatkan pendapatan yang lebih tinggi sebesar 0.97 kali dibandingkan umur sebelumnya. Hal ini menunjukkan semakin tua umur seseorang, maka kecenderungan untuk menghasilkan pendapatan cenderung menurun.
X6	Status Perkawinan	1.69	Responden dengan status menikah memiliki peluang 1.69 kali lebih besar dibanding responden dengan status belum menikah dalam mendapatkan pendapatan yang lebih besar.