regresi logistik Ordinal
Kelompok 8
2025-05-05
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
# Load library yang dibutuhkan
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── 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 errorslibrary(dplyr)
library(stringr)
library(ggplot2)
library(corrplot)## corrplot 0.95 loadedlibrary(smotefamily)
library(MASS)##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## selectlibrary(tibble)
library(tidyr)
library(rcompanion)
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBalibrary(caret)## Loading required package: lattice
##
## Attaching package: 'caret'
##
## The following object is masked from 'package:purrr':
##
## liftlibrary(DescTools)##
## Attaching package: 'DescTools'
##
## The following objects are masked from 'package:caret':
##
## MAE, RMSE# Baca dataset
data <- read.csv("survey.csv")
# Lihat 5 baris pertama
head(data)## Timestamp Age Gender Country state self_employed
## 1 2014-08-27 11:29:31 37 Female United States IL <NA>
## 2 2014-08-27 11:29:37 44 M United States IN <NA>
## 3 2014-08-27 11:29:44 32 Male Canada <NA> <NA>
## 4 2014-08-27 11:29:46 31 Male United Kingdom <NA> <NA>
## 5 2014-08-27 11:30:22 31 Male United States TX <NA>
## 6 2014-08-27 11:31:22 33 Male United States TN <NA>
## family_history treatment work_interfere no_employees remote_work
## 1 No Yes Often 6-25 No
## 2 No No Rarely More than 1000 No
## 3 No No Rarely 6-25 No
## 4 Yes Yes Often 26-100 No
## 5 No No Never 100-500 Yes
## 6 Yes No Sometimes 6-25 No
## tech_company benefits care_options wellness_program seek_help anonymity
## 1 Yes Yes Not sure No Yes Yes
## 2 No Don't know No Don't know Don't know Don't know
## 3 Yes No No No No Don't know
## 4 Yes No Yes No No No
## 5 Yes Yes No Don't know Don't know Don't know
## 6 Yes Yes Not sure No Don't know Don't know
## leave mental_health_consequence phys_health_consequence
## 1 Somewhat easy No No
## 2 Don't know Maybe No
## 3 Somewhat difficult No No
## 4 Somewhat difficult Yes Yes
## 5 Don't know No No
## 6 Don't know No No
## coworkers supervisor mental_health_interview phys_health_interview
## 1 Some of them Yes No Maybe
## 2 No No No No
## 3 Yes Yes Yes Yes
## 4 Some of them No Maybe Maybe
## 5 Some of them Yes Yes Yes
## 6 Yes Yes No Maybe
## mental_vs_physical obs_consequence comments
## 1 Yes No <NA>
## 2 Don't know No <NA>
## 3 No No <NA>
## 4 No Yes <NA>
## 5 Don't know No <NA>
## 6 Don't know No <NA>Tahap Preprocessing
Melihat beberapa baris awal (head(data)) untuk memahami struktur dan isi data secara kasat mata.
Melihat struktur data (str(data)) untuk untuk mengetahui tipe data setiap kolom, jumlah observasi (1259), dan variabel (27), serta membedakan mana yang bertipe numerik, karakter, atau faktor.
Melihat ringkasan statistik (summary(data)) untuk mendeteksi nilai ekstrim (outlier), distribusi nilai, dan nilai-nilai yang umum muncul.
Membersihkan kolom target work_interfere menghapus baris dengan nilai NA pada work_interfere karena kolom merupakan target prediksi dalam model
Membersihkan kolom usia (Age) menyimpan usia antara 18 hingga 100 tahun untuk memastikan validitas data demografis dan mencegah outlier memengaruhi model.
Transformasi Target Menjadi Faktor Ordinal Kolom work_interfere diubah menjadi faktor berurutan (ordinal factor): Never < Rarely < Sometimes < Often. karena akan digunakan dalam ordinal logistic regression, di mana urutan memiliki makna.
Menangani Missing Values Ditemukan beberapa NA, salah satunya pada kolom self_employed. karena memiliki banyak missing value serta tidak relevan untuk pemodelan, sehingga dihapus dari dataset.
Pembersihan dan Normalisasi Kolom Gender karena banyaknya variasi penulisan pada kolom gender sehingga dilakukan trim spasi dan mengkatagorikan ulang dengan case_when sehingga hanya terdapat “Male” dan “Female” pada kolom gender
Konversi Kolom Tipe Karakter Menjadi Faktor semua kolom yang masih bertipe karakter diubah menjadi tipe faktor, agar bisa digunakan dalam pemodelan. karena sebagian besar model di R mengharapkan data kategorik sebagai faktor, bukan karakter.
# Lihat struktur data
str(data)## 'data.frame': 1259 obs. of 27 variables:
## $ Timestamp : chr "2014-08-27 11:29:31" "2014-08-27 11:29:37" "2014-08-27 11:29:44" "2014-08-27 11:29:46" ...
## $ Age : num 37 44 32 31 31 33 35 39 42 23 ...
## $ Gender : chr "Female" "M" "Male" "Male" ...
## $ Country : chr "United States" "United States" "Canada" "United Kingdom" ...
## $ state : chr "IL" "IN" NA NA ...
## $ self_employed : chr NA NA NA NA ...
## $ family_history : chr "No" "No" "No" "Yes" ...
## $ treatment : chr "Yes" "No" "No" "Yes" ...
## $ work_interfere : chr "Often" "Rarely" "Rarely" "Often" ...
## $ no_employees : chr "6-25" "More than 1000" "6-25" "26-100" ...
## $ remote_work : chr "No" "No" "No" "No" ...
## $ tech_company : chr "Yes" "No" "Yes" "Yes" ...
## $ benefits : chr "Yes" "Don't know" "No" "No" ...
## $ care_options : chr "Not sure" "No" "No" "Yes" ...
## $ wellness_program : chr "No" "Don't know" "No" "No" ...
## $ seek_help : chr "Yes" "Don't know" "No" "No" ...
## $ anonymity : chr "Yes" "Don't know" "Don't know" "No" ...
## $ leave : chr "Somewhat easy" "Don't know" "Somewhat difficult" "Somewhat difficult" ...
## $ mental_health_consequence: chr "No" "Maybe" "No" "Yes" ...
## $ phys_health_consequence : chr "No" "No" "No" "Yes" ...
## $ coworkers : chr "Some of them" "No" "Yes" "Some of them" ...
## $ supervisor : chr "Yes" "No" "Yes" "No" ...
## $ mental_health_interview : chr "No" "No" "Yes" "Maybe" ...
## $ phys_health_interview : chr "Maybe" "No" "Yes" "Maybe" ...
## $ mental_vs_physical : chr "Yes" "Don't know" "No" "No" ...
## $ obs_consequence : chr "No" "No" "No" "Yes" ...
## $ comments : chr NA NA NA NA ...# Lihat ringkasan statistik
summary(data)## Timestamp Age Gender Country
## Length:1259 Min. :-1.726e+03 Length:1259 Length:1259
## Class :character 1st Qu.: 2.700e+01 Class :character Class :character
## Mode :character Median : 3.100e+01 Mode :character Mode :character
## Mean : 7.943e+07
## 3rd Qu.: 3.600e+01
## Max. : 1.000e+11
## state self_employed family_history treatment
## Length:1259 Length:1259 Length:1259 Length:1259
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## work_interfere no_employees remote_work tech_company
## Length:1259 Length:1259 Length:1259 Length:1259
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## benefits care_options wellness_program seek_help
## Length:1259 Length:1259 Length:1259 Length:1259
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## anonymity leave mental_health_consequence
## Length:1259 Length:1259 Length:1259
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## phys_health_consequence coworkers supervisor
## Length:1259 Length:1259 Length:1259
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## mental_health_interview phys_health_interview mental_vs_physical
## Length:1259 Length:1259 Length:1259
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## obs_consequence comments
## Length:1259 Length:1259
## Class :character Class :character
## Mode :character Mode :character
##
##
## # Lihat kolom target
table(data$work_interfere, useNA = "ifany")##
## Never Often Rarely Sometimes <NA>
## 213 144 173 465 264# Hapus baris dengan work_interfere NA
data_clean <- data %>% filter(!is.na(work_interfere))
# Cek nilai Age yang tidak masuk akal
summary(data_clean$Age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.726e+03 2.700e+01 3.100e+01 1.005e+08 3.600e+01 1.000e+11# Hapus age outlier (usia < 18 atau > 100 tidak masuk akal)
data_clean <- data_clean %>%
filter(Age >= 18 & Age <= 100)
# Konversi work_interfere menjadi faktor ordinal
data_clean$work_interfere <- factor(data_clean$work_interfere,
levels = c("Never", "Rarely", "Sometimes", "Often"),
ordered = TRUE)
# Cek hasil
table(data_clean$work_interfere)##
## Never Rarely Sometimes Often
## 212 173 464 140# Cek missing value per kolom
colSums(is.na(data_clean))## Timestamp Age Gender
## 0 0 0
## Country state self_employed
## 0 391 18
## family_history treatment work_interfere
## 0 0 0
## no_employees remote_work tech_company
## 0 0 0
## benefits care_options wellness_program
## 0 0 0
## seek_help anonymity leave
## 0 0 0
## mental_health_consequence phys_health_consequence coworkers
## 0 0 0
## supervisor mental_health_interview phys_health_interview
## 0 0 0
## mental_vs_physical obs_consequence comments
## 0 0 848# Lihat variasi Gender
table(tolower(trimws(data_clean$Gender))) # Trim dan lowercase untuk lihat bentuk aslinya##
## agender
## 1
## androgyne
## 1
## cis-female/femme
## 1
## cis female
## 1
## cis male
## 3
## cis man
## 1
## enby
## 1
## f
## 44
## femail
## 1
## female
## 161
## female (cis)
## 1
## female (trans)
## 2
## fluid
## 1
## genderqueer
## 1
## guy (-ish) ^_^
## 1
## m
## 110
## mail
## 1
## make
## 4
## mal
## 1
## male
## 638
## male-ish
## 1
## male leaning androgynous
## 1
## malr
## 1
## msle
## 1
## nah
## 1
## neuter
## 1
## non-binary
## 1
## ostensibly male, unsure what that really means
## 1
## queer
## 1
## queer/she/they
## 1
## trans-female
## 1
## trans woman
## 1
## woman
## 2# Normalisasi kolom Gender
data_clean$Gender <- tolower(trimws(data_clean$Gender)) # ubah ke huruf kecil dan buang spasi
data_clean$Gender <- case_when(
str_detect(data_clean$Gender, "female") ~ "Female",
str_detect(data_clean$Gender, "male") ~ "Male",
str_detect(data_clean$Gender, "trans") ~ "Trans",
TRUE ~ "Other"
)
# Cek hasil normalisasi Gender
table(data_clean$Gender)##
## Female Male Other Trans
## 167 644 177 1# Tangani missing value di self_employed
# karena proporsi NA kecil, diisi dengan "No" (karena sebagian besar orang tidak self-employed)
prop.table(table(data_clean$self_employed, useNA = "ifany"))##
## No Yes <NA>
## 0.8594540 0.1223458 0.0182002data_clean$self_employed[is.na(data_clean$self_employed)] <- "No"
library(dplyr)
# Drop kolom state dan comments karena tidak digunakan dalam model dan banyak NA
data_clean <- dplyr::select(data_clean, -state, -comments)
# Hanya simpan baris dengan Gender Male atau Female
data_clean <- data_clean %>%
filter(Gender %in% c("Male", "Female"))
# Cek jumlah akhir untuk memastikan
table(data_clean$Gender)##
## Female Male
## 167 644# Ubah semua kolom karakter jadi faktor, kecuali Age dan Timestamp
data_clean <- data_clean %>%
mutate(across(where(is.character) & !c(Timestamp), as.factor))EDA
Distribusi Target (work_interfere) distribusi kategori work_interfere dengan hasil menunjukkan jumlah observasi untuk masing-masing kategori: Never, Rarely, Sometimes, Often. visualisasi distribusi dengan bar chart membantu memahami ketidakseimbangan data antar kelas target, penting untuk pertimbangan SMOTE nanti. Target work_interfere diubah menjadi ordinal factor karena urutan kategori penting untuk pemodelan ordinal logistic regression.
Visualisasi Keterkaitan dengan Variabel Kategorikal semua variabel kategorikal (kecuali target) diubah ke format long agar bisa divisualisasikan proporsinya per variabel plot proporsi dibuat untuk tiap level work_interfere terhadap semua variabel kategorikal, yang bertujuan untuk melihat tren proporsi tiap level target pada masing-masing kategori prediktor
Ukur Hubungan antara Variabel Kategorikal dan Target menggunakan CramerV untuk mengukur kekuatan asosiasi antara work_interfere dan variabel karegorikal. Lalu hasil diurutkan dari yang paling tinggi, semakin tinggi nilainya (maks 1) akan semakin kuat aosiasinya
Visualisasi Variabel Numerik (Usia) terhadap Target boxplot digunakan untuk melihat distribusi usia (Age) terhadap setiap level work_interfere, untuk mengamati apakah ada perbedaan usia yang nyata antar level gangguan kerja (interferensi).
Uji Statistik: Hubungan antar Variabel Kategorikal dan Target setiap variabel kategorikal, dilakukan uji Chi-Square terhadap work_interfere. Lalu hasil akan berbentuk tabel: nama variabel, p-value, dan metode uji statistik yang digunakan.
Korelasi Antar Variabel Numerik (Setelah Label Encoding) variabel faktor dikonversi ke numerik (label encoding) secara manual agar bisa dihitung korelasinya.
``` r
#EDA
# Distribusi jumlah tiap kategori target
table(data_clean$work_interfere)##
## Never Rarely Sometimes Often
## 170 143 375 123# Visualisasi distribusi target
library(ggplot2)
ggplot(data_clean, aes(x = work_interfere)) +
geom_bar(fill = "#00BFC4") +
theme_minimal() +
labs(title = "Distribusi Target: work_interfere", x = "Work Interfere", y = "Jumlah")
# Ubah work_interfere menjadi ordered factor (urutan penting untuk regresi ordinal!)
data_clean$work_interfere <- factor(
data_clean$work_interfere,
levels = c("Never", "Rarely", "Sometimes", "Often"),
ordered = TRUE
)
# Ambil semua variabel faktor kecuali target
vars_cat_all <- names(data_clean)[sapply(data_clean, is.factor)]
vars_cat_all <- setdiff(vars_cat_all, "work_interfere")
# Data long format
library(tidyr)
library(dplyr)
data_long_all <- data_clean %>%
dplyr::select(all_of(vars_cat_all), work_interfere) %>%
pivot_longer(cols = all_of(vars_cat_all), names_to = "variable", values_to = "value")
# Buat folder untuk simpan gambar (opsional)
dir.create("eda_charts", showWarnings = FALSE)
# Loop untuk tiap variabel kategorikal
for (var_name in unique(data_long_all$variable)) {
# Filter hanya untuk 1 variabel
plot_data <- data_long_all[data_long_all$variable == var_name, ]
# Buat plot
p <- ggplot(plot_data, aes(x = value, fill = work_interfere)) +
geom_bar(position = "fill") +
theme_minimal() +
labs(
title = paste("Proporsi work_interfere terhadap", var_name),
x = NULL,
y = "Proporsi"
) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Tampilkan atau simpan
print(p) # tampilkan di viewer
# Simpan sebagai file (opsional)
ggsave(filename = paste0("eda_charts/", var_name, "_proporsi.png"), plot = p, width = 7, height = 5)
}
library(rcompanion)
# Hitung Cramér's V untuk semua variabel faktor vs target
cramer_results <- sapply(vars_cat_all, function(var) {
cramerV(table(data_clean[[var]], data_clean$work_interfere))
})
# Urutkan dari yang paling berpengaruh
sort(cramer_results, decreasing = TRUE)## treatment.Cramer V family_history.Cramer V
## 0.54990 0.28510
## Country.Cramer V obs_consequence.Cramer V
## 0.22170 0.15750
## care_options.Cramer V leave.Cramer V
## 0.12800 0.12080
## benefits.Cramer V seek_help.Cramer V
## 0.12050 0.11800
## mental_health_consequence.Cramer V Gender.Cramer V
## 0.11800 0.11580
## mental_vs_physical.Cramer V no_employees.Cramer V
## 0.10520 0.09798
## supervisor.Cramer V anonymity.Cramer V
## 0.09371 0.09169
## wellness_program.Cramer V remote_work.Cramer V
## 0.08821 0.08793
## mental_health_interview.Cramer V self_employed.Cramer V
## 0.08577 0.08351
## phys_health_consequence.Cramer V tech_company.Cramer V
## 0.06937 0.05756
## phys_health_interview.Cramer V coworkers.Cramer V
## 0.05738 0.05718# Boxplot
ggplot(data_clean, aes(x = work_interfere, y = Age, fill = work_interfere)) +
geom_boxplot() +
theme_minimal() +
labs(title = "Boxplot Usia terhadap work_interfere", x = "Work Interfere", y = "Age")
# Variabel terpilih dari hasil Cramér’s V dan boxplot
# Target: work_interfere (ordinal)
# Prediktor:
# - Kategorikal: treatment, family_history, Country, obs_consequence, care_options, leave
# - Numerik: Age
vars_used <- c("work_interfere", "treatment", "family_history", "Country",
"obs_consequence", "care_options", "leave", "Age")
# Subset data bersih
data_model <- data_clean[, vars_used]
# Pastikan variabel kategorikal difaktorkan
data_model$work_interfere <- factor(data_model$work_interfere, ordered = TRUE)
data_model$treatment <- as.factor(data_model$treatment)
data_model$family_history <- as.factor(data_model$family_history)
data_model$Country <- as.factor(data_model$Country)
data_model$obs_consequence <- as.factor(data_model$obs_consequence)
data_model$care_options <- as.factor(data_model$care_options)
data_model$leave <- as.factor(data_model$leave)
# Uji Chi-Square antar variabel kategorikal terhadap target
cat_vars <- c("treatment", "family_history", "Country", "obs_consequence", "care_options", "leave")
chi_results <- lapply(cat_vars, function(var) {
tbl <- table(data_model[[var]], data_model$work_interfere)
expected <- suppressWarnings(chisq.test(tbl)$expected)
if (any(expected < 5)) {
# Simulasikan p-value karena tabel besar
p <- fisher.test(tbl, simulate.p.value = TRUE, B = 10000)$p.value
method <- "Fisher (simulated)"
} else {
p <- chisq.test(tbl)$p.value
method <- "Chi-squared"
}
data.frame(Variable = var, p_value = p, Method = method)
}) %>% bind_rows()
print(chi_results)## Variable p_value Method
## 1 treatment 7.013191e-53 Chi-squared
## 2 family_history 3.219589e-14 Chi-squared
## 3 Country 3.072693e-01 Fisher (simulated)
## 4 obs_consequence 1.596437e-04 Chi-squared
## 5 care_options 1.728975e-04 Chi-squared
## 6 leave 3.910316e-04 Chi-squaredlibrary(smotefamily)
library(MASS)
library(tibble)
# Salin dataset
data_smote <- data_model
data_smote$work_interfere <- as.factor(data_smote$work_interfere)
# Pisahkan X dan y
X <- dplyr::select(data_smote, -work_interfere)
y <- data_smote$work_interfere
# Konversi semua faktor jadi numerik manual (label encoding)
X_encoded <- X %>%
mutate(across(where(is.factor), ~ as.numeric(as.factor(.))))
# Cek ulang apakah semua kolom numerik
str(X_encoded)## 'data.frame': 811 obs. of 7 variables:
## $ treatment : num 2 1 2 1 1 2 2 1 2 1 ...
## $ family_history : num 1 1 2 1 2 2 2 1 1 1 ...
## $ Country : num 37 7 36 37 37 37 37 7 37 6 ...
## $ obs_consequence: num 1 1 2 1 1 1 1 1 1 1 ...
## $ care_options : num 2 1 3 1 2 1 3 1 1 2 ...
## $ leave : num 3 2 2 1 1 2 4 1 1 1 ...
## $ Age : num 37 32 31 31 33 35 42 23 31 29 ...# Hitung matriks korelasi
cor_matrix <- cor(X_encoded)
# Install dan load package corrplot jika belum ada
# install.packages("corrplot")
library(corrplot)
# Visualisasi korelasi dengan pengaturan yang lebih baik
corrplot(cor_matrix,
method = "color", # Menggunakan warna untuk visualisasi
type = "upper", # Hanya menunjukkan bagian atas dari matriks (karena korelasi bersifat simetris)
tl.cex = 0.7, # Ukuran teks label fitur
tl.col = "black", # Warna teks label
addCoef.col = "black", # Menambahkan nilai korelasi di atas plot
number.cex = 0.7, # Ukuran teks angka korelasi
col = colorRampPalette(c("blue", "white", "red"))(200), # Palet warna dari biru ke merah
diag = FALSE) # Menyembunyikan diagonal (karena korelasi diagonal selalu 1)
Proses PCA
menggunakan fungsi (prcomp()) untuk memastikan setiap fitur distandarisasi (mean = 0, sd = 1) sebelum PCA. penting karena variabel hasil encoding mempunyai skala berbeda hasil PCA: - PC1 menjelaskan 22.25% dari total variasi. - PC1–PC3 secara kumulatif menjelaskan 54.25% variasi. - PC1–PC6 menjelaskan 90.57%, yang berarti 6 komponen cukup untuk mewakili hampir seluruh informasi dari data asli.
# PCA untuk reduksi dimensi
pca_result <- prcomp(X_encoded, scale. = TRUE)
summary(pca_result) # Lihat proporsi varians tiap komponen## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.2481 1.0667 1.0497 0.9955 0.9016 0.8596 0.81236
## Proportion of Variance 0.2225 0.1625 0.1574 0.1416 0.1161 0.1056 0.09428
## Cumulative Proportion 0.2225 0.3851 0.5425 0.6841 0.8002 0.9057 1.00000# Visualisasi PCA
library(factoextra)
fviz_eig(pca_result) # Scree plot
fviz_pca_var(pca_result, col.var = "contrib") # Kontribusi variabel
Analisis Modeling
1. Persiapan Data dan Pemodelan - Penanganan Imbalance Data: SMOTE diterapkan untuk menyeimbangkan distribusi kelas target (work_interfere). Never: 170
Rarely: 143
Sometimes: 375
Often: 369 - Pembagian Data: Data dibagi menjadi 80% training dan 20% testing dengan stratifikasi untuk mempertahankan distribusi kelas.
- Validasi Silang (5-fold CV): Model diuji dengan lima metode link function (cauchit, cloglog, logistic, loglog, probit). Hasilnya: Akurasi tertinggi: 0.4503 (metode cauchit), diikuti oleh logistic (0.4503) dan probit (0.4446).
Kappa tertinggi: 0.2078 (metode cauchit), menunjukkan kesepakatan model dengan data yang lemah (Kappa < 0.2 dianggap “poor”).
2. Hasil Model Regresi Logistik Ordinal (POLR) - Koefisien dan signifikansi statistik: Faktor Signifikan (p-value < 0.05):
Treatment (Mendapat perawatan kesehatan mental):
Koefisien: 1.895 (p < 0.001) → Efek positif kuat. Artinya, individu yang mendapat perawatan cenderung mengalami gangguan kesehatan mental yang lebih tinggi di tempat kerja.
Family History (Riwayat keluarga):
Koefisien: 0.379 (p = 0.0024) → Efek positif. Riwayat keluarga meningkatkan risiko gangguan.
Country (Negara asal):
Koefisien: -0.019 (p = 0.0004) → Efek negatif lemah. Negara dengan nilai lebih tinggi (misalnya, kode tertentu) mungkin terkait dengan gangguan lebih rendah.
Observed Consequence (Konsekuensi diskriminasi):
Koefisien: 0.427 (p = 0.0073) → Efek positif. Diskriminasi meningkatkan tingkat gangguan.
Leave (Kemudahan cuti kesehatan mental):
Koefisien: 0.123 (p = 0.0022) → Efek positif. Akses cuti yang lebih mudah dikaitkan dengan gangguan lebih tinggi (mungkin karena kebutuhan cuti yang lebih besar).
Faktor Tidak Signifikan:
Care Options (Opsi perawatan): Koefisien negatif (-0.092) tetapi tidak signifikan (p = 0.1875).
Age: Koefisien negatif (-0.006) dan tidak signifikan (p = 0.4667).
- Threshold (Intercepts):
Never|Rarely: 1.6939
Rarely|Sometimes: 2.6784
Sometimes|Often: 4.4240 Nilai threshold yang meningkat menunjukkan peningkatan kesulitan untuk beralih ke kategori gangguan yang lebih tinggi.
3. Evaluasi Performa Model
Akurasi: 48.1% menunjukkan model hanya sedikit lebih baik daripada prediksi acak berdasarkan kelas mayoritas.
Kappa: 0.2438
Sensitivitas dan Spesifisitas:
Kelas “Never”: Sensitivitas tinggi (73.5%) dan spesifisitas 90.3%, artinya model cukup baik dalam mengidentifikasi kasus tanpa gangguan.
Kelas “Rarely”: Sensitivitas 0% (model gagal memprediksi kategori ini).
Kelas “Sometimes” dan “Often”: Sensitivitas rendah (36% dan 67.1%), menunjukkan kesulitan model dalam membedakan kategori gangguan menengah hingga tinggi.
Mcnemar’s Test: P-value sangat kecil (1.315e-05)
4. Interpretasi Faktor Utama
Treatment dan Family History adalah prediktor terkuat.
Negara memiliki efek negatif lemah, mungkin terkait dengan perbedaan budaya atau kebijakan kesehatan mental di berbagai negara.
Diskriminasi (obs_consequence) dan kemudahan cuti (leave) memperburuk gangguan, mungkin karena stres tambahan atau kebutuhan cuti yang tidak terpenuhi.
5. Keterbatasan Model
Akurasi dan Kappa Rendah: Model memiliki kemampuan prediktif terbatas, karena variabel prediktor yang tidak mencukupi atau kurang relevan, kompleksitas hubungan faktor-faktor yang tidak tertangkap oleh model linier. Kelas “Rarely” Tidak Terprediksi: Indikasi ketidakmampuan model membedakan kategori gangguan rendah.
# Apply SMOTE
set.seed(123)
smote_result <- SMOTE(X_encoded, y, K = 5, dup_size = 2)
# Gabungkan hasil SMOTE
data_balanced <- smote_result$data
data_balanced$class <- factor(data_balanced$class,
levels = c("Never", "Rarely", "Sometimes", "Often"),
ordered = TRUE)
# Ganti nama target
colnames(data_balanced)[ncol(data_balanced)] <- "work_interfere"
# Cek hasil
table(data_balanced$work_interfere)##
## Never Rarely Sometimes Often
## 170 143 375 369set.seed(123)
trainIndex <- createDataPartition(data_balanced$work_interfere, p = 0.8, list = FALSE)
trainData <- data_balanced[trainIndex, ]
testData <- data_balanced[-trainIndex, ]
# Setup cross-validation
train_control <- trainControl(method = "cv", number = 5) # 5-fold CV
# Train model dengan polr dan CV
set.seed(123)
cv_model <- train(work_interfere ~ .,
data = data_balanced,
method = "polr",
trControl = train_control)
# Lihat hasil
print(cv_model)## Ordered Logistic or Probit Regression
##
## 1057 samples
## 7 predictor
## 4 classes: 'Never', 'Rarely', 'Sometimes', 'Often'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 845, 845, 846, 846, 846
## Resampling results across tuning parameters:
##
## method Accuracy Kappa
## cauchit 0.4503174 0.2078321
## cloglog 0.4171913 0.1184057
## logistic 0.4502906 0.2035866
## loglog 0.4399133 0.2037775
## probit 0.4446168 0.1935380
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was method = cauchit.#Model
model_polr <- polr(work_interfere ~ ., data = data_balanced, Hess = TRUE)
summary(model_polr)## Call:
## polr(formula = work_interfere ~ ., data = data_balanced, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## treatment 1.894604 0.145725 13.0012
## family_history 0.379008 0.124983 3.0325
## Country -0.019232 0.005446 -3.5311
## obs_consequence 0.427331 0.159182 2.6845
## care_options -0.092416 0.070112 -1.3181
## leave 0.123121 0.040219 3.0613
## Age -0.006114 0.008400 -0.7279
##
## Intercepts:
## Value Std. Error t value
## Never|Rarely 1.6939 0.4140 4.0920
## Rarely|Sometimes 2.6784 0.4210 6.3621
## Sometimes|Often 4.4240 0.4345 10.1830
##
## Residual Deviance: 2488.159
## AIC: 2508.159ctable <- coef(summary(model_polr))
# Hitung p-value manual:
p_values <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
# Gabungkan ke dalam tabel koefisien:
ctable <- cbind(ctable, "p value" = round(p_values, 4))
print(ctable)## Value Std. Error t value p value
## treatment 1.894604019 0.145725020 13.0012267 0.0000
## family_history 0.379007521 0.124983094 3.0324703 0.0024
## Country -0.019232111 0.005446420 -3.5311472 0.0004
## obs_consequence 0.427330775 0.159182377 2.6845357 0.0073
## care_options -0.092415728 0.070112278 -1.3181105 0.1875
## leave 0.123121334 0.040219055 3.0612687 0.0022
## Age -0.006114199 0.008400042 -0.7278771 0.4667
## Never|Rarely 1.693879837 0.413952354 4.0919681 0.0000
## Rarely|Sometimes 2.678412634 0.420992422 6.3621398 0.0000
## Sometimes|Often 4.424047729 0.434453477 10.1830183 0.0000predicted_classes <- predict(model_polr, newdata = testData, type = "class")
all_levels <- c("Never", "Rarely", "Sometimes", "Often")
predicted_classes <- factor(predicted_classes, levels = all_levels)
testData$work_interfere <- factor(testData$work_interfere, levels = all_levels)
# Lalu jalankan confusion matrix
conf_mat <- confusionMatrix(predicted_classes, testData$work_interfere)
print(conf_mat)## Confusion Matrix and Statistics
##
## Reference
## Prediction Never Rarely Sometimes Often
## Never 25 5 11 1
## Rarely 0 0 0 0
## Sometimes 7 13 27 23
## Often 2 10 37 49
##
## Overall Statistics
##
## Accuracy : 0.481
## 95% CI : (0.4117, 0.5508)
## No Information Rate : 0.3571
## P-Value [Acc > NIR] : 0.0001543
##
## Kappa : 0.2438
##
## Mcnemar's Test P-Value : 1.315e-05
##
## Statistics by Class:
##
## Class: Never Class: Rarely Class: Sometimes Class: Often
## Sensitivity 0.7353 0.0000 0.3600 0.6712
## Specificity 0.9034 1.0000 0.6815 0.6423
## Pos Pred Value 0.5952 NaN 0.3857 0.5000
## Neg Pred Value 0.9464 0.8667 0.6571 0.7857
## Prevalence 0.1619 0.1333 0.3571 0.3476
## Detection Rate 0.1190 0.0000 0.1286 0.2333
## Detection Prevalence 0.2000 0.0000 0.3333 0.4667
## Balanced Accuracy 0.8194 0.5000 0.5207 0.6568``` ## Including Plots
You can also embed plots, for example:
{r pressure, echo=FALSE} plot(pressure)
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.