1.1. Definisi

• Regresi logistik ordinal merupakan salah satu analisis regresi logitik yang digunakan untuk menganalisa hubungan antara variabel dependen dengan variabel independen, dimana variabel dependen bersifat polikotomus dengan skala ordinal.

• Polikotomus disini maksudnya adalah data kategorik di mana kategori nya berjumlah lebih dari dua.

1.2. Model Regresi Logistik Ordinal

Model yang dapat dipakai untuk regresi logistik ordinal adalah model logit.

Keterangan :

\(j = 1,2,...,j\) adalah kategori respon.

1.3. Langkah Analisis

Adapun langkah-langkah analisis regresi logistik ordinal adalah sebagai berikut :

1. Melakukan input data

2. Melakukan pengecekan asumsi multikolinieritas antar variabel bebas. Jika tidak terdapat multikolinieritas, maka analisis dapat dilanjutkan. Jika terdapat multikolinieritas, data terlebih dahulu harus dilakukan penanganan seperti transformasi data, penambahan data atau penghapusan pada salah satu variabel yang berkorelasi kuat dengan variabel lainnya.

3. Melakukan uji kecocokan model atau Goodness of Fit untuk mengetahui apakah model yang dibentuk sudah baik atau belum. Pengujian dapat dilakukan menggunakan Uji Ordinal Hosmer and Lemeshow, Pulksteni and Robinson, Pearson dan Deviance Chi Square, serta Uji Lipsitz. Jika diperoleh kesimpulan model cocok atau H0 diterima, maka analisis regresi logistik ordinal cocok untuk memodelkan data tersebut.

4. Setelah itu akan dilanjutkan dengan uji serentak menggunakan statistik uji G untuk mengetahui ada tidaknya pengaruh variabel prediktor terhadap variabel respon secara bersama-sama. Jika terdapat pengaruh atau H0 ditolak, maka lanjutkan dengan uji signifikansi parsial menggunaka uji Wald. Sebaliknya, jika tidak terdapat pengaruh, maka uji parsial tidak perlu dilakukan.

5. Terakhir lakukan interpretasi model menggunakan nilai Odds Ratio.

2.1. Deskripsi

Misal ingin diketahui apakah terdapat pengaruh antara umur, jenis kelamin serta status perkawinan terhadap tingkat partisipasi kerja bakti di suatu kelurahan. Tingkat partisipasi tersebut diukur dalam 3 kategori, yaitu Tidak Pernah, Jarang dan Sering. Karena tingkat partisipasi sebagai variabel dependen merupakan variabel kategorik berskala ordinal, maka regresi logistik ordinal merupakan metode yang sesuai untuk mencapai tujuan tersebut.

Berikut tabel deskripsi dari data yang digunakan.

2.2. Input data

• Langkah pertama yang harus dilakukan saat akan melakukan analisis data menggunakan R adalah melakukan input data.

• Banyak cara yang dapat dilakukan untuk mengimport data, disini saya menggunakan sintaks read_excel yang terdapat pada library (readxl)

library(readxl)

data <- read_excel("Data Latihan.xlsx", sheet = "Kerja Bakti")
data

## # A tibble: 199 × 4
##      age gender married kerjabakti
##    <dbl>  <dbl>   <dbl>      <dbl>
##  1    59      1       1          1
##  2    39      0       1          1
##  3    28      0       1          1
##  4    30      1       1          2
##  5    26      1       0          1
##  6    36      1       1          1
##  7    40      0       1          1
##  8    55      1       1          1
##  9    54      0       1          1
## 10    28      1       0          3
## # ℹ 189 more rows

Data tersebut terdiri dari 199 baris dan 4 kolom. Artinya, ada 199 responden yang akan digunakan sebagai sampel penelitian.

Sebelum dilakukan pemodelan, cek terlebih dahulu apakah skala data sudah sesuai dengan tabel deskripsi pada bagian awal. Ingat, terdapat 3 variabel kategorik pada data tersebut, yaitu tingkat partisipasi, jenis kelamin dan status perkawinan. Kamu dapat mengecek informasi tersebut melalui str()

#melihat sruktur data
str(data)

## tibble [199 × 4] (S3: tbl_df/tbl/data.frame)
##  $ age       : num [1:199] 59 39 28 30 26 36 40 55 54 28 ...
##  $ gender    : num [1:199] 1 0 0 1 1 1 0 1 0 1 ...
##  $ married   : num [1:199] 1 1 1 1 0 1 1 1 1 0 ...
##  $ kerjabakti: num [1:199] 1 1 1 2 1 1 1 1 1 3 ...

Dari output di atas, terlihat bahwa variabel age sudah benar karena merupakan data numerik, sedangkan ketiga variabel lain yang seharusnya berskala kategorik harus diubah terlebih dahulu menggunakan as.factor() agar menjadi data kategorik.

2.3. Mengubah data integer menjadi kategorik

#Mengubah data integer menjadi kategorik
data$gender<-as.factor(data$gender)
data$married<-as.factor(data$married)
data$kerjabakti<-as.factor(data$kerjabakti)

Kemudian periksa kembali skala data menggunakan str().

#memeriksa kembali sruktur data

str(data)

## tibble [199 × 4] (S3: tbl_df/tbl/data.frame)
##  $ age       : num [1:199] 59 39 28 30 26 36 40 55 54 28 ...
##  $ gender    : Factor w/ 2 levels "0","1": 2 1 1 2 2 2 1 2 1 2 ...
##  $ married   : Factor w/ 2 levels "0","1": 2 2 2 2 1 2 2 2 2 1 ...
##  $ kerjabakti: Factor w/ 3 levels "1","2","3": 1 1 1 2 1 1 1 1 1 3 ...

summary(data)

##       age       gender  married kerjabakti
##  Min.   :17.0   0:107   0: 60   1:74      
##  1st Qu.:25.0   1: 92   1:139   2:96      
##  Median :33.0                   3:29      
##  Mean   :35.3                             
##  3rd Qu.:44.0                             
##  Max.   :77.0

2.4. Eksplorasi

summary(data)

##       age       gender  married kerjabakti
##  Min.   :17.0   0:107   0: 60   1:74      
##  1st Qu.:25.0   1: 92   1:139   2:96      
##  Median :33.0                   3:29      
##  Mean   :35.3                             
##  3rd Qu.:44.0                             
##  Max.   :77.0

library(arsenal)
tab<-tableby(kerjabakti~ .,data = data)
summary(tab,text=TRUE)

## 
## 
## |             |    1 (N=74)     |    2 (N=96)     |    3 (N=29)     |  Total (N=199)  | p value|
## |:------------|:---------------:|:---------------:|:---------------:|:---------------:|-------:|
## |age          |                 |                 |                 |                 | < 0.001|
## |-  Mean (SD) | 40.811 (13.190) | 32.156 (12.056) | 31.621 (11.397) | 35.296 (13.053) |        |
## |-  Range     | 17.000 - 71.000 | 17.000 - 77.000 | 22.000 - 61.000 | 17.000 - 77.000 |        |
## |gender       |                 |                 |                 |                 |   0.809|
## |-  0         |   42 (56.8%)    |   50 (52.1%)    |   15 (51.7%)    |   107 (53.8%)   |        |
## |-  1         |   32 (43.2%)    |   46 (47.9%)    |   14 (48.3%)    |   92 (46.2%)    |        |
## |married      |                 |                 |                 |                 | < 0.001|
## |-  0         |   12 (16.2%)    |   30 (31.2%)    |   18 (62.1%)    |   60 (30.2%)    |        |
## |-  1         |   62 (83.8%)    |   66 (68.8%)    |   11 (37.9%)    |   139 (69.8%)   |        |

boxplot(data$age)

boxplot(age~kerjabakti, data)

2.5. Pembentukan model

Semua skala data sudah sesuai dengan tabel deskripsi di awal. Langkah selanjutnya adalah melakukan pemodelan regresi logistik ordinal menggunakan sintaks polr() yang terdapat pada library (MASS) .

#pembentukan model
library(MASS)
model_ord=polr(kerjabakti~age+gender+married, method='logistic',data=data)

summary(model_ord)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = kerjabakti ~ age + gender + married, data = data, 
##     method = "logistic")
## 
## Coefficients:
##             Value Std. Error t value
## age      -0.03629    0.01177  -3.084
## gender1   0.28602    0.27749   1.031
## married1 -1.04737    0.33355  -3.140
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -2.4827  0.4682    -5.3023
## 2|3  0.1174  0.4239     0.2769
## 
## Residual Deviance: 367.9638 
## AIC: 377.9638

• Untuk melakukan pemodelan ini, nama variabel dependen ditulis sebelum tanda ~ dan semua variabel independen ditulis setelah tanda ~. Method yang dipilih adalah logistik.

#uji asumsi multikolinieritas

library(car)

## Loading required package: carData

vif(model_ord)

## 
## Re-fitting to get Hessian

##      age   gender  married 
## 1.090172 1.006878 1.092370

#Uji GoF

library(generalhoslem)

## Loading required package: reshape

lipsitz.test(model_ord)

## 
##  Lipsitz goodness of fit test for ordinal response models
## 
## data:  formula:  kerjabakti ~ age + gender + married
## LR statistic = 13.096, df = 9, p-value = 0.1583

#Uji serentak
library(pscl)

## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis

pR2(model_ord)

## fitting null model for pseudo-r2

##           llh       llhNull            G2      McFadden          r2ML 
## -183.98190937 -199.03783777   30.11185680    0.07564355    0.14042385 
##          r2CU 
##    0.16239300

#Chisquare tabel
qchisq(0.95, 3)

## [1] 7.814728

#Uji Parsial
koef=coef(summary(model_ord))

## 
## Re-fitting to get Hessian

p_val_parsial=pnorm(abs(koef[,'t value']),lower.tail=FALSE)*2
tabel_uji_parsial=cbind(koef,'p value'=p_val_parsial)
tabel_uji_parsial

##                Value Std. Error    t value      p value
## age      -0.03629359 0.01176745 -3.0842346 2.040766e-03
## gender1   0.28601940 0.27749170  1.0307314 3.026668e-01
## married1 -1.04736526 0.33354681 -3.1400848 1.688989e-03
## 1|2      -2.48267904 0.46822257 -5.3023481 1.143226e-07
## 2|3       0.11737143 0.42394305  0.2768566 7.818902e-01

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

## 
## Re-fitting to get Hessian

#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
## age      -0.0363     0.0118 -3.0842 0.0020
## gender1   0.2860     0.2775  1.0307 0.3027
## married1 -1.0474     0.3335 -3.1401 0.0017
## 1|2      -2.4827     0.4682 -5.3023 0.0000
## 2|3       0.1174     0.4239  0.2769 0.7819

2.6. Pembentukan Model 2 (menghilangkan variabel gender)

model_ord2=polr(kerjabakti~age+married, method='logistic',data=data)

summary(model_ord2)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = kerjabakti ~ age + married, data = data, method = "logistic")
## 
## Coefficients:
##             Value Std. Error t value
## age      -0.03593    0.01178  -3.048
## married1 -1.02864    0.33247  -3.094
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -2.5860  0.4588    -5.6365
## 2|3  0.0023  0.4094     0.0057
## 
## Residual Deviance: 369.029 
## AIC: 377.029

2.7. Berikut model yang terbentuk.

1. Model untuk menghitung berapa peluang responden tidak pernah mengikuti kerja bakti

   

2. Model untuk menghitung berapa peluang responden tidak pernah dan jarang mengikuti kerja bakti

   

Kedua model tersebut nantinya dapat digunakan untuk memprediksi tingkat partisipasi seorang responden berdasarkan nilai peluang dari tiga kategori tersebut, sebagai simulasi akan dicontohkan sebagai berikut.

Misal seorang responden berusia 59 tahun, berjenis kelamin perempuan dan sudah menikah. Maka peluang masing-masing kategori untuk responden dengan kriteria tersebut dapat dihitung dengan cara berikut.

Jadi, peluang ia tidak pernah mengikuti kerja bakti sebesar 0.6034806.

Peluang ia jarang mengikuti kerja bakti sebesar 0.3499908.

Sedangkan, peluang ia sering mengikuti kerja bakti adalah sebesar 0.0465286.

Karena dari ketiga kategori yang paling besar adalah peluang tidak pernah mengikuti kerja bakti, maka responden dengan karakteristik berusia 59 tahun, berjenis kelamin perempuan dan sudah menikah diprediksi tidak pernah mengikuti kerja bakti.

Jika Menggunakan syntax di R, hasilnya akan sama

#Input data baru (usia 59, perempuan, sudah menikah)
new = data.frame(age = 59, 
                 gender = as.factor("1"), 
                 married = as.factor("1")
                 )

prediksi <- predict(model_ord,newdata = new,"probs")
peluang <- data.frame(cumulatif = cumsum(prediksi),point=prediksi) 
peluang

##   cumulatif      point
## 1 0.6034805 0.60348051
## 2 0.9534713 0.34999075
## 3 1.0000000 0.04652874

2.8. Odds ratio

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

#Odds Ratio
exp(coef(model_ord))

##       age   gender1  married1 
## 0.9643571 1.3311183 0.3508610

Dari output tersebut, diperoleh informasi sebagai berikut.

1. Untuk variabel umur, didapatkan nilai odds ratio sebesar 0.96, artinya jika umur responden bertambah 1 tahun, maka kecenderungan ia untuk semakin terlibat kerja bakti sebesar 0.96 kali dibandingkan umur sebelumnya. Hal ini menunjukkan semakin tua umur seseorang, maka kecenderungan untuk terlibat akan semakin menurun.

2. Untuk variabel jenis kelamin, didapatkan nilai odds ratio sebesar 1.33, artinya responden perempuan memiliki kecenderungan 1.33 kali lebih besar untuk semakin terlibat dalam kerja bakti dibandingkan responden pria.

3. Untuk variabel status perkawinan, didapatkan nilai odds ratio sebesar 0.35, artinya responden yang sudah menikah memiliki kecenderungan 0.35 kali untuk semakin terlibat dalam kerja bakti dibandingkan responden yang belum menikah.