Regresi Logistik Ordinal Contoh 2

1. Package

library(ggplot2)
library(tidyverse)
library(gmodels)
library(ggmosaic)
library(corrplot)
library(caret)
library(rpart)
library(rpart.plot)
library(C50)
library(gmodels)
library(tidyr)
library(dplyr)
library(knitr)
library(rsample)
library(readxl)
library(MASS)

2. Input Data

data_akreditasi<-read_xlsx("Data Latihan.xlsx", sheet = "Akreditasi")
data_akreditasi

# A tibble: 104 × 4
     UN1   UN2   UN3 Akreditasi
   <dbl> <dbl> <dbl>      <dbl>
 1  68.1  68.6  71.8          4
 2  72.7  75.7  30.2          4
 3  46.4  48.5  46.4          4
 4  66.1  67.3  74.9          3
 5  54.3  64.4  75.4          4
 6  48.5  42    36.9          3
 7  70.5  76.3  40.6          4
 8  70.9  58.5  63.6          4
 9  59.7  38.2  65.9          3
10  50.6  68.5  53.0          3
# ℹ 94 more rows

#melihat sruktur data
str(data_akreditasi)

tibble [104 × 4] (S3: tbl_df/tbl/data.frame)
 $ UN1       : num [1:104] 68.1 72.7 46.4 66.1 54.3 48.5 70.5 70.9 59.7 50.6 ...
 $ UN2       : num [1:104] 68.6 75.7 48.5 67.3 64.4 ...
 $ UN3       : num [1:104] 71.8 30.2 46.4 74.9 75.4 ...
 $ Akreditasi: num [1:104] 4 4 4 3 4 3 4 4 3 3 ...

2.1. Mengubah data integer menjadi kategorik

#Mengubah data integer menjadi kategorik
data_akreditasi$Akreditasi<-as.factor(data_akreditasi$Akreditasi)

data_akreditasi

# A tibble: 104 × 4
     UN1   UN2   UN3 Akreditasi
   <dbl> <dbl> <dbl> <fct>     
 1  68.1  68.6  71.8 4         
 2  72.7  75.7  30.2 4         
 3  46.4  48.5  46.4 4         
 4  66.1  67.3  74.9 3         
 5  54.3  64.4  75.4 4         
 6  48.5  42    36.9 3         
 7  70.5  76.3  40.6 4         
 8  70.9  58.5  63.6 4         
 9  59.7  38.2  65.9 3         
10  50.6  68.5  53.0 3         
# ℹ 94 more rows

str(data_akreditasi)

tibble [104 × 4] (S3: tbl_df/tbl/data.frame)
 $ UN1       : num [1:104] 68.1 72.7 46.4 66.1 54.3 48.5 70.5 70.9 59.7 50.6 ...
 $ UN2       : num [1:104] 68.6 75.7 48.5 67.3 64.4 ...
 $ UN3       : num [1:104] 71.8 30.2 46.4 74.9 75.4 ...
 $ Akreditasi: Factor w/ 4 levels "1","2","3","4": 4 4 4 3 4 3 4 4 3 3 ...

summary(data_akreditasi)

      UN1             UN2             UN3        Akreditasi
 Min.   :36.10   Min.   :29.00   Min.   :27.20   1:16      
 1st Qu.:47.45   1st Qu.:44.53   1st Qu.:44.02   2: 3      
 Median :60.73   Median :52.20   Median :52.81   3:21      
 Mean   :58.85   Mean   :55.86   Mean   :54.41   4:64      
 3rd Qu.:68.03   3rd Qu.:68.33   3rd Qu.:65.09             
 Max.   :83.39   Max.   :86.45   Max.   :82.04

3. Pembentukan model

library(MASS)
model_akreditasi=polr(Akreditasi~., method='logistic',data=data_akreditasi)

3.1. Uji asumsi multikolinieritas

#uji asumsi multikolinieritas

library(car)
vif(model_akreditasi)

     UN1      UN2      UN3 
2.058265 1.949834 1.314553

3.2. Uji kesesuaian model.

#Uji GoF

library(generalhoslem)
lipsitz.test(model_akreditasi)


    Lipsitz goodness of fit test for ordinal response models

data:  formula:  Akreditasi ~ UN1 + UN2 + UN3
LR statistic = 14.137, df = 9, p-value = 0.1175

3.3. Uji serentak

#Uji serentak
library(pscl)
pR2(model_akreditasi)

fitting null model for pseudo-r2

          llh       llhNull            G2      McFadden          r2ML 
-100.93016529 -105.25590856    8.65148654    0.04109739    0.07982128 
         r2CU 
   0.09197122

#Chisquare tabel
qchisq(0.95, 2)

[1] 5.991465

library(lmtest) 
lrtest(model_akreditasi)

Likelihood ratio test

Model 1: Akreditasi ~ UN1 + UN2 + UN3
Model 2: Akreditasi ~ 1
  #Df  LogLik Df  Chisq Pr(>Chisq)  
1   6 -100.93                       
2   3 -105.26 -3 8.6515     0.0343 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.4. Uji Parsial

#Uji Parsial
koef=coef(summary(model_akreditasi))
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
UN1 -0.01438293 0.02419240 -0.5945227 0.55216258
UN2  0.04390920 0.02068839  2.1224078 0.03380351
UN3  0.01186448 0.01749116  0.6783130 0.49757326
1|2  0.46206163 1.08862203  0.4244463 0.67124033
2|3  0.68722838 1.09118411  0.6298006 0.52882507
3|4  1.78213263 1.09729112  1.6241202 0.10435017

4. Pembentukan Model 2

model_akreditasi2=polr(Akreditasi~UN2, method='logistic',data=data_akreditasi)

4.1. Uji Parsial

#Uji Parsial
koef2=coef(summary(model_akreditasi2))
p_val_parsial2=pnorm(abs(koef2[,'t value']),lower.tail=FALSE)*2
tabel_uji_parsial2=cbind(koef2,'p value'=p_val_parsial2)
tabel_uji_parsial2

         Value Std. Error   t value     p value
UN2 0.04043132 0.01488468 2.7163047 0.006601514
1|2 0.47992661 0.82620793 0.5808787 0.561322195
2|3 0.70402595 0.82727662 0.8510164 0.394760258
3|4 1.79162648 0.83915931 2.1350254 0.032758932

summary(model_akreditasi)

Call:
polr(formula = Akreditasi ~ ., data = data_akreditasi, method = "logistic")

Coefficients:
       Value Std. Error t value
UN1 -0.01438    0.02419 -0.5945
UN2  0.04391    0.02069  2.1224
UN3  0.01186    0.01749  0.6783

Intercepts:
    Value   Std. Error t value
1|2  0.4621  1.0886     0.4244
2|3  0.6872  1.0912     0.6298
3|4  1.7821  1.0973     1.6241

Residual Deviance: 201.8603 
AIC: 213.8603

summary(model_akreditasi2)

Call:
polr(formula = Akreditasi ~ UN2, data = data_akreditasi, method = "logistic")

Coefficients:
      Value Std. Error t value
UN2 0.04043    0.01488   2.716

Intercepts:
    Value  Std. Error t value
1|2 0.4799 0.8262     0.5809 
2|3 0.7040 0.8273     0.8510 
3|4 1.7916 0.8392     2.1350 

Residual Deviance: 202.4944 
AIC: 210.4944

5. Odds ratio

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

#Odds Ratio
exp(coef(model_akreditasi2))

    UN2 
1.04126

Jika nilai UN2 tahunsebelumnya mengalami kenaikan 1 point, kecenderungan ia untuk mendapatkan nilai akreditasi yang lebih tinggi sebesar 1.04 kali dibandingkan nilai sebelumnya. Hal ini menunjukkan semakin tinggi nilai UN2, maka kecenderungan untuk mendapatkan nilai akreditasi yang lebih baik semakin tinggi.