Analisis Regresi Logistik Ordinal untuk FeedbackScore Pelanggan

Dataset yang kami gunakan berisi tentang kepuasan customer berdasarkan demografi yang bermacam macam dan faktor kebiasaan. Dataset ini kami peroleh dari kaggle (link tercantum dibawah) dengan kolom kolom seperti CustomerID (berisi ID untuk customer), Age (umur), Gender (jenis kelamin), Country (negara), Income (pendapatan), ProductQuality (kualitas dari produk), ServiceQuality (kualitas pelayanan), PurchaseFrequency (frekuensi pembelian), FeedbackScore (skor kepuasan pelanggan), LoyaltyLevel (tingkat kesetiaan pelanggan), SatisfactionSCore (skor kepuasan pelanggan). Dalam dataset ini, terdapat variabel dependen berupa FeedBackScore yang merupakan variabel bertipe ordinal (Low, Medium, dan High). Tujuan kami melakukan analisis ini adalah ingin melakukan analisis guna mengetahui pengaruh faktor-faktor(Age, Income, ProductQuality, ServiceQuality, PurchaseFrequency, LoyaltyLevel,dan SatisfactionScore) terhadap FeedBackScore.

Link Kaggle : https://www.kaggle.com/datasets/jahnavipaliwal/customer-feedback-and-satisfaction/data

#Load Dataset

data <- read.csv("C:/Kuliah/semester 4/Analisis Multivariat/tugas/data/Dataset logistik ordinal.csv")
 
head(data)
##   CustomerID Age Gender Country Income ProductQuality ServiceQuality
## 1          1  56   Male      UK  83094              5              8
## 2          2  69   Male      UK  86860             10              2
## 3          3  46 Female     USA  60173              8             10
## 4          4  32 Female      UK  73884              7             10
## 5          5  60   Male      UK  97546              6              4
## 6          6  25   Male  France  49364              5              7
##   PurchaseFrequency FeedbackScore LoyaltyLevel SatisfactionScore
## 1                 5           Low       Bronze            100.00
## 2                 8        Medium         Gold            100.00
## 3                18        Medium       Silver            100.00
## 4                16           Low         Gold            100.00
## 5                13           Low       Bronze             82.00
## 6                19          High       Silver             80.71
str(data)
## 'data.frame':    38444 obs. of  11 variables:
##  $ CustomerID       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Age              : int  56 69 46 32 60 25 38 56 36 40 ...
##  $ Gender           : chr  "Male" "Male" "Female" "Female" ...
##  $ Country          : chr  "UK" "UK" "USA" "UK" ...
##  $ Income           : int  83094 86860 60173 73884 97546 49364 39716 45305 98001 85307 ...
##  $ ProductQuality   : int  5 10 8 7 6 5 10 10 10 6 ...
##  $ ServiceQuality   : int  8 2 10 10 4 7 6 10 2 4 ...
##  $ PurchaseFrequency: int  5 8 18 16 13 19 11 10 12 20 ...
##  $ FeedbackScore    : chr  "Low" "Medium" "Medium" "Low" ...
##  $ LoyaltyLevel     : chr  "Bronze" "Gold" "Silver" "Gold" ...
##  $ SatisfactionScore: num  100 100 100 100 82 ...

#Data Preprocessing

sapply(data, function(x) any(is.infinite(x)))
##        CustomerID               Age            Gender           Country 
##             FALSE             FALSE             FALSE             FALSE 
##            Income    ProductQuality    ServiceQuality PurchaseFrequency 
##             FALSE             FALSE             FALSE             FALSE 
##     FeedbackScore      LoyaltyLevel SatisfactionScore 
##             FALSE             FALSE             FALSE

Cek data infinite

colSums(is.na(data[, c("Age", "Income", "ProductQuality", "ServiceQuality", 
                       "PurchaseFrequency", "LoyaltyLevel", "SatisfactionScore")]))
##               Age            Income    ProductQuality    ServiceQuality 
##                 0                 0                 0                 0 
## PurchaseFrequency      LoyaltyLevel SatisfactionScore 
##                 0                 0                 0

Cek data NA

data$FeedbackScore <- factor(data$FeedbackScore, 
                             levels = c("Low", "Medium", "High"), 
                             ordered = TRUE)

Pastikan FeedBackScore telah ordinal

#Explore Data

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.2
ggplot(data, aes(x = FeedbackScore)) +
  geom_bar(fill = "skyblue") +
  theme_minimal()

Cek persebaran data pada variabel dependen (FeedBackScore)
Dari barchart di atas terlihat persebaran kategori balance, dan karena balance, dapat langsung dilakukan modelling tanpa perlu penyesuaian khusus seperti undersampling atau oversampling.

ggplot(data, aes(x = Age)) +
  geom_histogram(binwidth = 5, fill = "lightblue", color = "black") +
  theme_minimal() +
  labs(title = "Distribusi Usia Pelanggan", x = "Age", y = "Jumlah")

Cek distribusi usia pelanggan

Terlihat bahwa jumlah pelanggan dalam rentang usia 20-65 tahun terlihat sampir sama rata, artinya populasi pelanggan tersebar cukup merata di berbagai kelompok umur.

ggplot(data, aes(x = FeedbackScore, y = Income, fill = FeedbackScore)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Income berdasarkan FeedbackScore", x = "FeedbackScore", y = "Income")

Cek income berdasarkan FeedBackScore

Terlihat ketiga feedback(low, medium, high) memiliki median yang hampir sejajar di sekitar 75,000. Rentang IQR pun juga hampir sama, ini menunjukkan bahwa variasi penghasilan antar pelanggan tiap kategori hampir seragam. Maka dari itu, dapat disimpulkan bahwa income bukan faktor dominan dalam FeedBackScore

Uji Validitas dan Reliabilitas

# Hitung R Tabel
n <- nrow(data)
df <- n - 2
r_table <- qt(1-0.05/2, df) / sqrt(qt(1-0.05/2, df)^2 + df)
print(r_table)
## [1] 0.009996249
# Variabel yang diuji validitasnya
variabels <- c("Age", "Income", "ProductQuality", "ServiceQuality", "PurchaseFrequency")

# Inisialisasi list variabel valid
variabel_valid <- c()

# Loop untuk uji korelasi masing-masing variabel terhadap SatisfactionScore
for (var in variabels) {
  cat("\nKorelasi antara", var, "dengan SatisfactionScore:\n")
  hasil <- cor.test(data[[var]], data$SatisfactionScore)
  r_hitung <- abs(as.numeric(hasil$estimate))
  print(paste("r hitung:", r_hitung))
  print(paste("p-value:", hasil$p.value))
  
  # Cek validitas
  if (r_hitung >= r_table) {
    variabel_valid <- c(variabel_valid, var)
  }
}
## 
## Korelasi antara Age dengan SatisfactionScore:
## [1] "r hitung: 0.157510374539876"
## [1] "p-value: 5.01388016192675e-212"
## 
## Korelasi antara Income dengan SatisfactionScore:
## [1] "r hitung: 0.24512907004616"
## [1] "p-value: 0"
## 
## Korelasi antara ProductQuality dengan SatisfactionScore:
## [1] "r hitung: 0.547690440857993"
## [1] "p-value: 0"
## 
## Korelasi antara ServiceQuality dengan SatisfactionScore:
## [1] "r hitung: 0.553614228283434"
## [1] "p-value: 0"
## 
## Korelasi antara PurchaseFrequency dengan SatisfactionScore:
## [1] "r hitung: 0.113018173028152"
## [1] "p-value: 1.75483783206418e-109"
# Tampilkan variabel yang valid
cat("\nVariabel yang valid berdasarkan r tabel:\n")
## 
## Variabel yang valid berdasarkan r tabel:
print(variabel_valid)
## [1] "Age"               "Income"            "ProductQuality"   
## [4] "ServiceQuality"    "PurchaseFrequency"
# Uji Reliabilitas

library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
cek <- alpha(data[, c("Age", "Income","ProductQuality", "ServiceQuality", "PurchaseFrequency", "SatisfactionScore")])

cek
## 
## Reliability analysis   
## Call: alpha(x = data[, c("Age", "Income", "ProductQuality", "ServiceQuality", 
##     "PurchaseFrequency", "SatisfactionScore")])
## 
##   raw_alpha std.alpha G6(smc) average_r  S/N     ase  mean   sd median_r
##    0.00038      0.42    0.57      0.11 0.73 1.3e-05 12538 4330   0.0051
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt    -0.02     0  0.02
## Duhachek  0.00     0  0.00
## 
##  Reliability if an item is dropped:
##                   raw_alpha std.alpha G6(smc) average_r    S/N alpha se   var.r
## Age                 4.0e-04    0.4639  0.6002   0.14756 0.8655  1.0e-05 5.1e-02
## Income              3.4e-01    0.4436  0.5679   0.13750 0.7971  4.1e-03 5.1e-02
## ProductQuality      4.0e-04    0.3777  0.4258   0.10827 0.6071  1.3e-05 3.2e-02
## ServiceQuality      4.0e-04    0.3697  0.4210   0.10499 0.5865  1.3e-05 3.2e-02
## PurchaseFrequency   4.0e-04    0.4703  0.6138   0.15082 0.8880  1.3e-05 5.2e-02
## SatisfactionScore   1.5e-06    0.0039  0.0032   0.00079 0.0039  8.1e-06 2.4e-05
##                     med.r
## Age               0.00566
## Income            0.00566
## ProductQuality    0.00562
## ServiceQuality    0.00064
## PurchaseFrequency 0.00517
## SatisfactionScore 0.00064
## 
##  Item statistics 
##                       n    raw.r std.r r.cor   r.drop    mean      sd
## Age               38444  0.00058  0.38 0.100  4.3e-06    43.5    15.0
## Income            38444  1.00000  0.41 0.181  1.5e-01 75076.6 25975.8
## ProductQuality    38444 -0.00115  0.51 0.455 -1.3e-03     5.5     2.9
## ServiceQuality    38444  0.00560  0.52 0.470  5.5e-03     5.5     2.9
## PurchaseFrequency 38444  0.00151  0.37 0.068  1.3e-03    10.5     5.8
## SatisfactionScore 38444  0.24598  0.86 1.013  2.5e-01    85.3    16.9
## 
## Non missing response frequency for each item
##                  1   2   3   4   5   6   7   8   9  10 miss
## ProductQuality 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1    0
## ServiceQuality 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1    0

Berdasarkan hasil validitas dan reliabilitas tersebut didapatkan bahwasanya data tersebut teruji secara validitas karena nilai r ≥ r tabel, tetapi tidak teruji secara reliabilitas karena nilai Cronbach’s Alpha < 0.6.

#Uji Multikolinearitas

library(car)  # untuk fungsi vif()
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:psych':
## 
##     logit
# Buat formula untuk regresi
formula_vif <- as.formula(
  paste("SatisfactionScore ~", paste(variabel_valid, collapse = " + "))
)

# Fit model regresi linier
model_vif <- lm(formula_vif, data = data)

# Hitung VIF
vif_values <- vif(model_vif)

print(vif_values)
##               Age            Income    ProductQuality    ServiceQuality 
##          1.000124          1.000030          1.000115          1.000115 
## PurchaseFrequency 
##          1.000053

Berdasarkan hasil dari uji multikolinearitas dengan perhitungan VIF.Didapatkan semua variabel independent dikatakan bebas dari multikolinearitas (VIF < 10).

#Distribusi Ordinal Pada FeedBackScore

library(ggplot2)

# Hitung jumlah masing-masing kategori
feedback_counts <- table(data$FeedbackScore)

# Ubah ke data frame supaya mudah dipakai ggplot
feedback_df <- as.data.frame(feedback_counts)
colnames(feedback_df) <- c("FeedbackScore", "Count")

# Buat pie chart
ggplot(feedback_df, aes(x = "", y = Count, fill = FeedbackScore)) +
  geom_col(width = 1) +
  coord_polar(theta = "y") + 
  theme_void() + 
  labs(title = "Distribusi Ordinal FeedbackScore") +
  scale_fill_brewer(palette = "Set3")

Berdasarkan visualisasi distribusi ordinal pada FeedBackScore terlihat bahwa jumlahnya sama, sehingga persebaran data untuk setiap ordinal itu seimbang

#Pemodelan

library(MASS)
model <- polr(FeedbackScore ~Age + ProductQuality + ServiceQuality + 
                PurchaseFrequency + LoyaltyLevel + SatisfactionScore, 
              data = data, Hess=TRUE) # Hapus Kolom Income karena mengandung data infinite
summary(model)
## Call:
## polr(formula = FeedbackScore ~ Age + ProductQuality + ServiceQuality + 
##     PurchaseFrequency + LoyaltyLevel + SatisfactionScore, data = data, 
##     Hess = TRUE)
## 
## Coefficients:
##                         Value Std. Error  t value
## Age                -2.042e-04  0.0006473 -0.31556
## ProductQuality     -4.101e-03  0.0044110 -0.92961
## ServiceQuality      2.484e-03  0.0044087  0.56352
## PurchaseFrequency   3.523e-05  0.0016517  0.02133
## LoyaltyLevelGold   -4.660e-03  0.0229171 -0.20334
## LoyaltyLevelSilver -5.491e-03  0.0229962 -0.23878
## SatisfactionScore  -7.151e-04  0.0009235 -0.77427
## 
## Intercepts:
##             Value    Std. Error t value 
## Low|Medium   -0.7785   0.0570   -13.6543
## Medium|High   0.5994   0.0570    10.5239
## 
## Residual Deviance: 84463.98 
## AIC: 84481.98

Didapatkan model dengan Interpretasi Koefisien :
Berdasarkan koefisienya didapatkan insight bahwasanya :
- Jika koefisien positif, maka semakin besar nilai variabel itu, semakin besar juga kemungkinan responden memilih kategori yang lebih tinggi (kategori Medium atau High).
- Jika koefisien negatif, maka semakin besar nilai variabel itu, semakin besar juga kemungkinan responden memilih kategori yang lebih rendah (kategori Low atau Medium).

#Uji G ( Likelihood Ratio Test)

library(MASS)

# Model full
model_full <- polr(FeedbackScore ~ Age + ProductQuality + ServiceQuality + 
                     PurchaseFrequency + LoyaltyLevel + SatisfactionScore, 
                   data = data, Hess = TRUE)

# Model null
model_null <- polr(FeedbackScore ~ 1, data = data, Hess = TRUE)

# Hitung G
G_statistic <- 2 * (logLik(model_full) - logLik(model_null))
G_statistic
## 'log Lik.' 4.804971 (df=9)
# Hitung derajat bebas
df_g <- length(coef(model_full))

# Hitung p-value
p_value <- pchisq(G_statistic, df = df_g, lower.tail = FALSE)

# Hitung chi-square tabel (x^2 tabel)
chi_sq_table <- qchisq(0.95, df = df_g)

# Output
cat("G Statistic (hitung):", G_statistic, "\n")
## G Statistic (hitung): 4.804971
cat("Degrees of Freedom:", df_g, "\n")
## Degrees of Freedom: 7
cat("Chi-Square Tabel (5%):", chi_sq_table, "\n")
## Chi-Square Tabel (5%): 14.06714
cat("P-Value:", p_value, "\n")
## P-Value: 0.6837495

Berdasarkan hasil uji G didapakan bahwasanya nilai Statistik G lebih kecil daripada nilai tabel chi square. Sehingga gagal tolak H0 dimana H0 itu menandakan bahwasanya model yang dibuat tidak ada variabel independen yang mempengaruhi variabel dependen yakni FeedbackScore.

#Kesimpulan

Dari hasil pengujian validitas dan reliabilitas instrumen penelitian, diperoleh hasi bahwa uji validitas terpenuhi namun uji reliabilitas tidak terpenuhi. Lalu variabel indenpende diharapkan diuji Multikolineritas, sehingga didapatkan hasil bahwa antar variabel independent tidak terdapat korelasi. Dengan kata lain, tidak terjadi kasus multikolinearitas.
Model logistik ordinal yang dibangun menggunakan metode Proportional Odds Model (polr) menunjukkan bahwa prediktor yang digunakan tidak mampu secara signifikan mempengaruhi tingkat FeedbackScore.
Meskipun uji validitas sudah menunjukkan bahwa variabel-variabel independen layak diukur, namun ketidakterpenuhinya reliabilitas dan hasil uji G yang tidak signifikan mengindikasikan bahwa model tersebut belum memadai untuk menjelaskan variasi dalam tingkat kepuasan pengguna.

Dengan demikian, disarankan untuk:

Refrensi :
1. Winarko, M. T. D., & Kartini, A. Y. . (2022). Analisis Kepuasan Pengguna Jasa Petugas Parkir Dinas Perhubungan Bojonegoro Menggunakan Regresi Logistik Ordinal. Jurnal Statistika Dan Komputasi, 1(1), 31–41. https://doi.org/10.32665/statkom.v1i1.442
2. Hidayah, N., Indahwati, Fitrianto, A., Erfiani, & Aliu, M. A. (2024). Pemodelan tingkat kecanduan games online menggunakan regresi logistik ordinal. MATH LOCUS: Jurnal Riset dan Inovasi Pendidikan Matematika, 5(1). https://doi.org/10.31002/mathlocus.v5i1.4335