Persiapan

Memuat paket

library(MASS)     # untuk analisis diskriminan: fungsi lda() dan qda()

## Warning: package 'MASS' was built under R version 4.5.3

library(dplyr)    # untuk manipulasi data: group_by(), summarise(), dll.

## Warning: package 'dplyr' was built under R version 4.5.3

## 
## Attaching package: 'dplyr'

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

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)  # untuk membuat grafik

## Warning: package 'ggplot2' was built under R version 4.5.3

library(knitr)    # untuk menampilkan tabel rapi dengan kable()

Fungsi bantu dan warna

Beberapa fungsi kecil dibuat di awal agar kode berikutnya singkat dan konsisten.

id() membulatkan angka dan menampilkannya dengan koma desimal.
fmt_p() menampilkan p-value; nilai yang sangat kecil ditulis < 0,0001.
sd_pooled() menghitung simpangan baku gabungan dalam kelompok untuk membakukan koefisien diskriminan.
kable_d() menampilkan tabel dengan jumlah desimal tetap pada kolom tertentu (misalnya 0,0190, bukan 0,019).
tabel_konfusi() membuat matriks klasifikasi (aktual vs prediksi) dalam bentuk tabel rapi.
warna menetapkan satu warna tetap untuk tiap status agar semua grafik seragam.

# Format angka dibulatkan 'd' digit, desimal berupa koma
id <- function(x, d = 2) format(round(x, d), nsmall = d, decimal.mark = ",")

# Format p-value: p yang sangat kecil ditulis "< 0,0001"
fmt_p <- function(p) ifelse(p < 0.0001, "< 0,0001", id(p, 4))

# Simpangan baku gabungan (pooled within-group) tiap peubah
sd_pooled <- function(X, g) {
  W <- Reduce(`+`, lapply(levels(g), function(l) {
    Xl <- X[g == l, , drop = FALSE]
    (nrow(Xl) - 1) * cov(Xl)                 
  }))
  sqrt(diag(W / (nrow(X) - nlevels(g))))     
}

# Matriks klasifikasi: baris = status aktual, kolom = status prediksi
tabel_konfusi <- function(aktual, prediksi) {
  m   <- table(Aktual = aktual, Prediksi = prediksi)
  out <- data.frame(rownames(m), unclass(m), check.names = FALSE, row.names = NULL)
  names(out)[1] <- "Aktual \u2193 / Prediksi \u2192"
  out
}

# Tabel dengan kolom tertentu diformat memakai jumlah desimal tetap (koma desimal)
kable_d <- function(df, kolom, d = 4, ...) {
  df[kolom] <- lapply(df[kolom], id, d = d)   
  kable(df, ...)
}

# Warna tetap untuk setiap status
warna <- c("Lanjut" = "#1b9e77", "Drop Out" = "#d95f02", "Pindah" = "#7570b3")

Impor Data

Sumber data

Data diambil dari tabel pada slide Studi Kasus, yang memuat 85 mahasiswa. Tabel pada slide tersusun dalam beberapa blok kolom, sehingga data dimasukkan per status dengan urutan baca dari kiri ke kanan, atas ke bawah:

Tiga data frame (lanjut, drop_out, pindah) dibuat dengan data.frame(). Setiap data frame berisi vektor IPK dan Skor (dibuat dengan c(...)) serta kolom Status.
rbind() menumpuk ketiganya menjadi satu data dat.
factor(..., levels = ...) mengubah Status menjadi faktor dengan urutan tetap (Lanjut, Drop Out, Pindah).
X (matriks peubah), g (kelompok), p (jumlah peubah), dan k (jumlah kelompok) disiapkan untuk dipakai di langkah selanjutnya.

# Status "Lanjut" (31 mahasiswa) 
lanjut <- data.frame(
  IPK  = c(2.96, 3.14, 3.22, 3.29, 3.69, 3.46, 3.03, 3.19, 3.63, 3.59, 3.3, 3.4,
           3.5, 3.78, 3.44, 3.48, 3.47, 3.35, 3.39, 3.28, 3.21, 3.58,
           3.33, 3.4, 3.38, 3.26, 3.6, 3.37, 3.8, 3.76, 3.24),
  Skor = c(596, 473, 482, 527, 505, 693, 626, 663, 447, 588, 563, 553,
           572, 591, 692, 528, 552, 520, 543, 523, 530, 564,
           565, 431, 605, 664, 609, 559, 521, 646, 467),
  Status = "Lanjut")

# Status "Drop Out" (28 mahasiswa) 
drop_out <- data.frame(
  IPK  = c(2.54, 2.43, 2.2, 2.36, 2.57, 2.35, 2.51, 2.51, 2.36, 2.36, 2.66, 2.68, 2.48,
           2.46, 2.63, 2.44, 2.13, 2.41, 2.55, 2.31, 2.41, 2.19, 2.35, 2.6, 2.55,
           2.72, 2.85, 2.9),
  Skor = c(446, 425, 474, 531, 542, 406, 412, 458, 399, 482, 420, 414, 533,
           509, 504, 336, 408, 469, 538, 505, 489, 411, 321, 394, 528,
           399, 381, 384),
  Status = "Drop Out")

# Status "Pindah" (26 mahasiswa) 
pindah <- data.frame(
  IPK  = c(2.86, 2.85, 3.14, 3.28, 2.89, 3.15, 3.5,
           2.89, 2.8, 3.13, 3.01, 2.79, 2.89, 2.91, 2.75, 2.73, 3.12, 3.08,
           3.03, 3, 3.03, 3.05, 2.85, 3.01, 3.03, 3.04),
  Skor = c(494, 496, 419, 371, 447, 313, 402,
           485, 444, 416, 471, 490, 431, 446, 546, 467, 463, 440,
           419, 509, 438, 399, 483, 453, 414, 446),
  Status = "Pindah")

dat <- rbind(lanjut, drop_out, pindah)                       
dat$Status <- factor(dat$Status, levels = c("Lanjut", "Drop Out", "Pindah"))  

X <- as.matrix(dat[, c("IPK", "Skor")])   
g <- dat$Status                           
p <- ncol(X)                              
k <- nlevels(g)                           
n_kelompok <- table(g)                    

# Pemeriksaan kelengkapan data: harus 85 baris dan tanpa data hilang
stopifnot(nrow(dat) == 85, !anyNA(dat))

write.csv(dat, "data_studi_kasus_s2.csv", row.names = FALSE)

kable(head(dat, 6), caption = "Enam baris pertama data")

Enam baris pertama data
IPK	Skor	Status
2.96	596	Lanjut
3.14	473	Lanjut
3.22	482	Lanjut
3.29	527	Lanjut
3.69	505	Lanjut
3.46	693	Lanjut

kable(as.data.frame(n_kelompok), col.names = c("Status", "Jumlah"),
      caption = "Jumlah mahasiswa per status")

Jumlah mahasiswa per status
Status	Jumlah
Lanjut	31
Drop Out	28
Pindah	26

Data berisi 85 mahasiswa: 31 berstatus Lanjut, 28 Drop Out, dan 26 Pindah. Jumlah tiap kelompok relatif seimbang (26 sampai 31), jadi tidak ada kelompok yang terlalu mendominasi. Tidak ada data hilang (pemeriksaan stopifnot() lolos). IPK berada pada rentang 2,13 sampai 3,80 dan skor tes pada 313 sampai 693, nilai yang wajar untuk skala IPK 0 sampai 4.

Eksplorasi Data

Ringkasan per status

group_by(Status) memisahkan data menurut status, lalu summarise() menghitung rata-rata, simpangan baku, minimum, dan maksimum untuk kedua peubah. Prefiks across(..., list(...)) menghasilkan semua statistik sekaligus.

ringkas <- dat %>%
  group_by(Status) %>%                                 
  summarise(
    `Rata-rata IPK`  = mean(IPK),
    `SB IPK`         = sd(IPK),
    `Min IPK`        = min(IPK),
    `Maks IPK`       = max(IPK),
    `Rata-rata Skor` = mean(Skor),
    `SB Skor`        = sd(Skor),
    `Min Skor`       = min(Skor),
    `Maks Skor`      = max(Skor),
    .groups = "drop"
  )
kable(ringkas, digits = c(0, 3, 3, 2, 2, 1, 1, 0, 0),
      caption = "Ringkasan IPK dan skor tes per status (SB = simpangan baku)")

Ringkasan IPK dan skor tes per status (SB = simpangan baku)
Status	Rata-rata IPK	SB IPK	Min IPK	Maks IPK	Rata-rata Skor	SB Skor	Min Skor	Maks Skor
Lanjut	3.404	0.209	2.96	3.8	561.2	68.0	431	693
Drop Out	2.482	0.183	2.13	2.9	447.1	62.4	321	542
Pindah	2.993	0.172	2.73	3.5	446.2	47.4	313	546

Lanjut: IPK rata-rata 3,40 dan skor tes rata-rata 561,2; keduanya tertinggi.
Drop Out: IPK rata-rata 2,48 (terendah) dan skor tes 447,1.
Pindah: IPK rata-rata 2,99 berada di antara Lanjut dan Drop Out, tetapi skor tes rata-ratanya (446,2) hampir sama dengan Drop Out.

Pola ini memberi petunjuk awal: skor tes cenderung memisahkan Lanjut dari dua kelompok lain, sedangkan IPK memisahkan ketiganya secara bertingkat.

Sebaran tiap peubah

Data diubah ke bentuk panjang dengan rbind() agar IPK dan Skor dapat digambar berdampingan dalam dua panel (facet_wrap()) dengan skala sumbu bebas (scales = "free_y") karena satuannya berbeda. geom_boxplot() menggambar kotak sebaran, sedangkan geom_jitter(height = 0) menambahkan titik amatan tanpa menggeser nilai vertikalnya.

# Membuat data format panjang
panjang <- rbind(
  data.frame(Status = dat$Status, Peubah = "IPK S1", Nilai = dat$IPK),
  data.frame(Status = dat$Status, Peubah = "Skor tes masuk S2", Nilai = dat$Skor)
)

# Menyamakan penulisan kategori
levels(panjang$Status)[levels(panjang$Status) == "Drop out"] <- "Drop Out"

# Membuat boxplot
ggplot(panjang, aes(x = Status, y = Nilai, fill = Status)) +
  geom_boxplot(alpha = 0.65, outlier.shape = NA, show.legend = FALSE) +
  geom_jitter(width = 0.12, height = 0, size = 1.4, alpha = 0.55, show.legend = FALSE) +
  scale_fill_manual(values = warna) +
  facet_wrap(~ Peubah, scales = "free_y") +
  labs(title = "Sebaran peubah menurut status", x = NULL, y = "Nilai") +
  theme_minimal(base_size = 11)

Pada IPK, ketiga kotak tersusun bertingkat (Lanjut paling tinggi, Drop Out paling rendah, Pindah di tengah) dengan sedikit tumpang tindih. Pada skor tes, kotak Lanjut berada jelas di atas, sedangkan kotak Drop Out dan Pindah hampir berimpit. Artinya skor tes sulit membedakan Drop Out dari Pindah, tetapi cukup baik membedakan Lanjut dari keduanya.

Hubungan IPK dan skor tes

Diagram pencar menempatkan IPK pada sumbu-x dan skor tes pada sumbu-y; warna dan bentuk titik menunjukkan status. stat_ellipse() menggambar elips yang melingkupi sekitar 90% sebaran tiap status, sehingga bentuk dan ukuran sebaran masing-masing kelompok mudah dibandingkan. cor() menghitung korelasi IPK dan skor tes di dalam tiap status.

ggplot(dat, aes(x = IPK, y = Skor, colour = Status, shape = Status)) +
  geom_point(size = 2.8, alpha = 0.85) +
  stat_ellipse(level = 0.90, linewidth = 0.7) +             
  scale_colour_manual(values = warna) +
  labs(title = "Sebaran IPK S1 dan skor tes masuk S2 menurut status",
       x = "IPK S1", y = "Skor tes masuk S2") +
  theme_minimal(base_size = 12)

# Korelasi IPK dan skor tes di dalam tiap status
korelasi <- sapply(levels(g), function(l) cor(X[g == l, "IPK"], X[g == l, "Skor"]))
kable(data.frame(Status = names(korelasi), `Korelasi IPK-Skor` = round(korelasi, 3),
                 check.names = FALSE, row.names = NULL),
      caption = "Korelasi IPK dan skor tes di dalam tiap status")

Korelasi IPK dan skor tes di dalam tiap status
Status	Korelasi IPK-Skor
Lanjut	0.004
Drop Out	-0.104
Pindah	-0.662

Kelompok Lanjut (hijau) berada di kanan-atas dan terpisah cukup jelas, terutama pada sumbu skor tes. Drop Out (oranye) berada di kiri-bawah. Pindah (ungu) berada di tengah pada sumbu IPK dan tumpang tindih dengan Drop Out pada sumbu skor tes. Elips ketiga kelompok memiliki bentuk dan kemiringan yang tidak seragam, yang mengisyaratkan matriks ragam-peragam yang berbeda (akan diuji formal berikutnya). Korelasi IPK-skor di dalam kelompok berbeda-beda (Lanjut 0,00, Drop Out -0,10, Pindah -0,66) dan paling kuat pada kelompok Pindah. Pola hubungan kedua peubah yang berlainan antarkelompok ini kembali mengisyaratkan matriks ragam-peragam yang tidak sama.

Uji Asumsi

Dua asumsi utama analisis diskriminan diperiksa pada seluruh data (85 mahasiswa): normalitas multivariat di tiap kelompok dan kesamaan matriks ragam-peragam antarkelompok. Satu pemeriksaan tambahan, yaitu pencilan multivariat, membantu menjelaskan hasilnya.

Uji normalitas multivariat (Mardia)

Uji Mardia memeriksa dua hal pada tiap kelompok:

Skewness (kemencengan) multivariat. Statistik b1 diuji dengan sebaran khi-kuadrat.
Kurtosis (keruncingan) multivariat. Statistik b2 diubah menjadi nilai z dan diuji dengan sebaran normal baku.

Di dalam fungsi mardia(): scale(Y, scale = FALSE) memusatkan data pada rata-ratanya, D adalah matriks jarak Mahalanobis antar-amatan, lalu b1 dan b2 dihitung dari D. Fungsi dijalankan untuk tiap status dengan lapply() dan hasilnya digabung dengan rbind().

Hipotesis: H₀ = data kelompok tersebut menyebar normal ganda; H₁ = tidak normal ganda. Taraf nyata 5%.

mardia <- function(Y) {
  Y  <- as.matrix(Y)
  n  <- nrow(Y);  p <- ncol(Y)
  Yc <- scale(Y, scale = FALSE)                 
  S  <- cov(Y) * (n - 1) / n                    
  D  <- Yc %*% solve(S) %*% t(Yc)               
  b1 <- sum(D^3) / n^2                          
  b2 <- mean(diag(D)^2)                         
  z_kurt <- (b2 - p * (p + 2)) / sqrt(8 * p * (p + 2) / n)  
  data.frame(
    n = n,
    b1 = b1,
    p_skewness = pchisq(n * b1 / 6, df = p * (p + 1) * (p + 2) / 6, lower.tail = FALSE),
    b2 = b2,
    z_kurtosis = z_kurt,
    p_kurtosis = 2 * pnorm(-abs(z_kurt))
  )
}

# Jalankan untuk tiap status lalu gabungkan
hasil_mardia <- do.call(rbind, lapply(levels(g), function(l) {
  cbind(Status = l, mardia(X[g == l, , drop = FALSE]))
}))

# Keputusan: tolak H0 bila salah satu p-value < 0,05
hasil_mardia$Keputusan <- ifelse(pmin(hasil_mardia$p_skewness, hasil_mardia$p_kurtosis) < 0.05,
                                 "Tolak H0", "Gagal tolak H0")

kable(hasil_mardia, digits = c(0, 0, 3, 4, 3, 3, 4, 0),
      col.names = c("Status", "n", "Skewness (b1)", "p Skewness", "Kurtosis (b2)",
                    "z Kurtosis", "p Kurtosis", "Keputusan"),
      caption = "Uji normalitas multivariat Mardia per status")

Uji normalitas multivariat Mardia per status
Status	n	Skewness (b1)	p Skewness	Kurtosis (b2)	z Kurtosis	p Kurtosis	Keputusan
Lanjut	31	0.091	0.9762	6.827	-0.816	0.4144	Gagal tolak H0
Drop Out	28	0.815	0.4330	6.515	-0.982	0.3260	Gagal tolak H0
Pindah	26	1.855	0.0901	11.188	2.032	0.0422	Tolak H0

H₀ tidak ditolak pada 2 kelompok, tetapi ditolak pada kelompok Pindah (ada p-value di bawah 0,05), sehingga asumsi normal ganda tidak terpenuhi pada kelompok tersebut. p-value terkecilnya adalah 0,0422, jadi pelanggarannya tergolong ringan (mendekati 0,05) dan kelompok tersebut hanya berisi 26 amatan, sehingga hasil klasifikasi untuk kelompok itu perlu dibaca lebih hati-hati. Analisis diskriminan umumnya cukup tahan terhadap pelanggaran ringan asumsi normalitas, terutama bila ukuran kelompok tidak terlalu kecil.

Pemeriksaan pencilan multivariat

Pencilan multivariat adalah mahasiswa yang kombinasi IPK dan skornya jauh dari pusat kelompoknya. mahalanobis() menghitung jarak kuadrat Mahalanobis (D²) tiap mahasiswa terhadap rata-rata kelompoknya, dengan memperhitungkan ragam dan korelasi kedua peubah. Nilai D² dibandingkan dengan titik kritis khi-kuadrat qchisq(0.975, df = 2); mahasiswa dengan D² di atas titik tersebut ditandai.

jarak <- do.call(rbind, lapply(levels(g), function(l) {
  Xl <- X[g == l, , drop = FALSE]
  data.frame(Status = l,
             Baris  = which(g == l),                              
             IPK    = Xl[, "IPK"], Skor = Xl[, "Skor"],
             D2     = mahalanobis(Xl, colMeans(Xl), cov(Xl)))      
}))
p <- ncol(X)
batas_d2 <- qchisq(0.975, df = p)                                  
pencilan <- subset(jarak, D2 > batas_d2)                           

cat("Titik kritis D2 (khi-kuadrat 97,5%, db = 2):", id(batas_d2, 2), "\n")

## Titik kritis D2 (khi-kuadrat 97,5%, db = 2): 7,38

if (nrow(pencilan) > 0) {
  kable(pencilan, digits = c(0, 0, 2, 0, 2), row.names = FALSE,
        caption = "Mahasiswa dengan jarak Mahalanobis melebihi titik kritis")
} else cat("Tidak ada pencilan multivariat.\n")

Mahasiswa dengan jarak Mahalanobis melebihi titik kritis
Status	Baris	IPK	Skor	D2
Pindah	65	3.15	313	9.49
Pindah	66	3.50	402	10.50

Terdapat 2 mahasiswa yang melewati titik kritis (D² > 7,38), semuanya berada pada kelompok Pindah. Mereka adalah kombinasi IPK-skor yang tidak lazim bagi kelompoknya, misalnya skor tes sangat rendah dibanding IPK-nya. Sebagai pemeriksaan, p-value kurtosis kelompok Pindah berubah dari 0,0422 menjadi 0,5655 setelah pencilan tersebut dikeluarkan. Jadi pelanggaran normalitas pada kelompok ini sebagian besar disebabkan oleh beberapa pencilan, bukan oleh bentuk sebaran secara keseluruhan. Data tidak dibuang pada analisis utama karena tidak ada alasan substantif bahwa pencatatannya keliru.

Uji kesamaan matriks ragam-peragam (Box’s M)

Uji Box’s M membandingkan matriks ragam-peragam antarkelompok.

Si adalah daftar matriks kovarians tiap kelompok, Sp adalah matriks kovarians gabungan (pooled).
M membandingkan determinan matriks gabungan dengan determinan tiap kelompok; semakin besar M, semakin berbeda matriks-matriksnya.
C adalah faktor koreksi agar chi = C × M mendekati sebaran khi-kuadrat dengan derajat bebas df = p(p+1)(k-1)/2.

Hipotesis: H₀ = matriks ragam-peragam semua kelompok sama; H₁ = minimal satu berbeda. Taraf nyata 5%.

box_m <- function(X, g) {
  gs <- levels(g)                  
  n  <- table(g)                   
  N  <- sum(n);  k <- length(gs);  p <- ncol(X)
  Si <- lapply(gs, function(l) cov(X[g == l, , drop = FALSE]))                        
  Sp <- Reduce(`+`, Map(function(S, ni) (ni - 1) * S, Si, as.list(n))) / (N - k)       
  M  <- (N - k) * log(det(Sp)) - sum((n - 1) * sapply(Si, function(S) log(det(S))))
  u  <- sum(1 / (n - 1)) - 1 / (N - k)
  C  <- 1 - (2 * p^2 + 3 * p - 1) / (6 * (p + 1) * (k - 1)) * u                      
  chi <- C * M
  df  <- p * (p + 1) * (k - 1) / 2
  list(M = M, chi = chi, df = df, p = pchisq(chi, df, lower.tail = FALSE), Si = setNames(Si, gs))
}

bm <- box_m(X, g)

kable(data.frame(
  Statistik = c("M", "Khi-kuadrat terkoreksi", "Derajat bebas", "p-value"),
  Nilai     = c(id(bm$M, 3), id(bm$chi, 3), bm$df, fmt_p(bm$p))
), caption = "Uji Box's M untuk kesamaan matriks ragam-peragam")

Uji Box’s M untuk kesamaan matriks ragam-peragam
Statistik	Nilai
M	16,665
Khi-kuadrat terkoreksi	16,074
Derajat bebas	6
p-value	0,0134

# Ragam tiap peubah dan kovarians IPK-skor per kelompok (komponen matriks kovarians)
kable(data.frame(
  Status = names(bm$Si),
  `Ragam IPK`      = sapply(bm$Si, function(S) S["IPK", "IPK"]),
  `Ragam Skor`     = sapply(bm$Si, function(S) S["Skor", "Skor"]),
  `Kovarians IPK-Skor` = sapply(bm$Si, function(S) S["IPK", "Skor"]),
  check.names = FALSE, row.names = NULL),
  digits = 3, caption = "Komponen matriks ragam-peragam tiap kelompok")

Komponen matriks ragam-peragam tiap kelompok
Status	Ragam IPK	Ragam Skor	Kovarians IPK-Skor
Lanjut	0.044	4618.247	0.058
Drop Out	0.034	3891.254	-1.192
Pindah	0.030	2246.905	-5.404

Statistik Box’s M terkoreksi bernilai 16,074 dengan derajat bebas 6 dan p-value 0,0134. Karena p-value < 0,05, H₀ ditolak: matriks ragam-peragam ketiga kelompok tidak sama. Tabel komponen memperlihatkan bahwa ragam skor dan kovarians IPK-skor memang berbeda antarkelompok. Akibatnya asumsi dasar LDA tidak terpenuhi dan QDA lebih sesuai secara teori; LDA tetap dibangun sebagai pembanding.

Box’s M peka terhadap pelanggaran normalitas dan cenderung mudah menolak H₀. Karena itu pilihan LDA atau QDA tidak hanya ditentukan oleh uji ini, melainkan juga dibuktikan lewat kinerja prediksi pada bagian Evaluasi.

Ringkasan uji asumsi

kable(data.frame(
  Asumsi = c("Normal ganda (Mardia)", "Pencilan multivariat (Mahalanobis)",
             "Kesamaan matriks ragam-peragam (Box's M)"),
  Hasil = c(ifelse(n_tolak_mardia == 0, "Tidak ada kelompok yang menolak H0",
                   paste0("H0 ditolak pada: ", paste(grup_tolak, collapse = ", "))),
            paste0(nrow(pencilan), " mahasiswa melewati titik kritis"),
            paste0("p ", fmt_p(bm$p))),
  Kesimpulan = c(ifelse(n_tolak_mardia == 0, "Terpenuhi", "Terpenuhi sebagian (pelanggaran ringan)"),
                 ifelse(nrow(pencilan) == 0, "Tidak ada", "Perlu dicermati"),
                 ifelse(bm$p < 0.05, "Tidak terpenuhi (utamakan QDA)", "Terpenuhi (LDA sesuai)"))
), caption = "Ringkasan pemeriksaan asumsi")

Ringkasan pemeriksaan asumsi
Asumsi	Hasil	Kesimpulan
Normal ganda (Mardia)	H0 ditolak pada: Pindah	Terpenuhi sebagian (pelanggaran ringan)
Pencilan multivariat (Mahalanobis)	2 mahasiswa melewati titik kritis	Perlu dicermati
Kesamaan matriks ragam-peragam (Box’s M)	p 0,0134	Tidak terpenuhi (utamakan QDA)

Uji Perbedaan Antarkelompok

Sebelum membangun fungsi diskriminan, perlu dipastikan bahwa ketiga status memang berbeda pada kedua peubah.

manova(cbind(IPK, Skor) ~ Status) menguji kedua peubah sekaligus. Statistik Wilks’ Lambda mendekati 0 bila kelompok sangat terpisah dan mendekati 1 bila tidak terpisah.
aov() (ANOVA satu arah) menguji tiap peubah sendiri-sendiri.
TukeyHSD() adalah uji lanjut yang membandingkan pasangan status, agar terlihat pasangan mana yang berbeda.

# MANOVA: uji serentak kedua peubah
fit_manova <- manova(cbind(IPK, Skor) ~ Status, data = dat)
st <- summary(fit_manova, test = "Wilks")$stats
kable(data.frame(
  Sumber          = c("Status", "Galat (Residuals)"),
  `Derajat bebas` = c(st[1, "Df"], st[2, "Df"]),
  `Wilks' Lambda` = c(id(st[1, "Wilks"], 4), ""),
  `F hampiran`    = c(id(st[1, "approx F"], 3), ""),
  `df1; df2`      = c(paste0(st[1, "num Df"], "; ", st[1, "den Df"]), ""),
  `p-value`       = c(fmt_p(st[1, "Pr(>F)"]), ""),
  check.names = FALSE),
  caption = "MANOVA (uji Wilks' Lambda)")

MANOVA (uji Wilks’ Lambda)
Sumber	Derajat bebas	Wilks’ Lambda	F hampiran	df1; df2	p-value
Status	2	0,1264	73,426	4; 162	< 0,0001
Galat (Residuals)	82

# ANOVA satu arah untuk tiap peubah
anova_ipk  <- aov(IPK  ~ Status, data = dat)
anova_skor <- aov(Skor ~ Status, data = dat)

tabel_anova <- function(fit) {
  a <- summary(fit)[[1]]
  data.frame(Sumber = c("Status", "Galat"), `Derajat bebas` = a[["Df"]],
             `Jumlah kuadrat` = a[["Sum Sq"]], `Kuadrat tengah` = a[["Mean Sq"]],
             F = c(a[["F value"]][1], NA), `p-value` = c(fmt_p(a[["Pr(>F)"]][1]), NA),
             check.names = FALSE)
}
kable(tabel_anova(anova_ipk),  digits = 3, caption = "ANOVA untuk IPK S1 menurut status")

ANOVA untuk IPK S1 menurut status
Sumber	Derajat bebas	Jumlah kuadrat	Kuadrat tengah	F	p-value
Status	2	12.502	6.251	173.305	< 0,0001
Galat	82	2.958	0.036	NA	NA

kable(tabel_anova(anova_skor), digits = 1, caption = "ANOVA untuk skor tes masuk S2 menurut status")

ANOVA untuk skor tes masuk S2 menurut status
Sumber	Derajat bebas	Jumlah kuadrat	Kuadrat tengah	F	p-value
Status	2	258471.1	129235.6	35.3	< 0,0001
Galat	82	299783.9	3655.9	NA	NA

# Uji lanjut Tukey 
tabel_tukey <- function(fit) {
  t <- TukeyHSD(fit)$Status
  data.frame(Pasangan = rownames(t), Selisih = t[, "diff"], `Batas bawah` = t[, "lwr"],
             `Batas atas` = t[, "upr"], `p (disesuaikan)` = fmt_p(t[, "p adj"]),
             check.names = FALSE, row.names = NULL)
}
kable(tabel_tukey(anova_ipk),  digits = 3, caption = "Uji lanjut Tukey untuk IPK S1")

Uji lanjut Tukey untuk IPK S1
Pasangan	Selisih	Batas bawah	Batas atas	p (disesuaikan)
Drop Out-Lanjut	-0.921	-1.040	-0.803	< 0,0001
Pindah-Lanjut	-0.411	-0.532	-0.291	< 0,0001
Pindah-Drop Out	0.510	0.387	0.634	< 0,0001

kable(tabel_tukey(anova_skor), digits = 1, caption = "Uji lanjut Tukey untuk skor tes masuk S2")

Uji lanjut Tukey untuk skor tes masuk S2
Pasangan	Selisih	Batas bawah	Batas atas	p (disesuaikan)
Drop Out-Lanjut	-114.2	-151.8	-76.5	< 0,0001
Pindah-Lanjut	-115.0	-153.4	-76.6	< 0,0001
Pindah-Drop Out	-0.8	-40.1	38.5	0,9986

Wilks’ Lambda = 0,1264 (p < 0,0001). Nilai yang jauh di bawah 1 dan p yang sangat kecil berarti vektor rata-rata ketiga status berbeda nyata, sehingga analisis diskriminan layak dilakukan.

ANOVA: IPK (F = 173,3) dan skor tes (F = 35,3) sama-sama berbeda nyata antarstatus. Nilai F IPK jauh lebih besar, jadi IPK adalah pembeda keseluruhan yang lebih kuat.
Tukey (IPK): ketiga pasangan status berbeda nyata, sehingga IPK membedakan ketiga status secara bertingkat (Lanjut > Pindah > Drop Out).
Tukey (skor tes): pasangan yang tidak berbeda nyata: Pindah-Drop Out (p = 0,9986). Skor tes hanya memisahkan Lanjut dari dua status lain, tetapi tidak dapat membedakan Drop Out dari Pindah.

Kesimpulan sementara: IPK memisahkan ketiga status, sedangkan skor tes terutama memisahkan Lanjut dari yang tidak lanjut. Kombinasi keduanya diharapkan menghasilkan klasifikasi terbaik.

Pembangunan Model Diskriminan

Model LDA

lda(Status ~ IPK + Skor, data = dat) membangun fungsi diskriminan linear dengan Status sebagai kelompok dan IPK serta Skor sebagai prediktor.

Karena ada 3 kelompok, terbentuk 2 fungsi diskriminan (jumlah kelompok − 1).
Koefisien asli sulit dibandingkan karena satuan IPK (0-4) dan skor (ratusan) berbeda jauh. Karena itu koefisien dibakukan dengan mengalikannya dengan simpangan baku gabungan (sd_pooled()).
predict(lda_fit)$x menghasilkan skor LD1 dan LD2 setiap mahasiswa; rata-ratanya per status disebut sentroid kelompok.

lda_fit <- lda(Status ~ IPK + Skor, data = dat)   
lda_fit

## Call:
## lda(Status ~ IPK + Skor, data = dat)
## 
## Prior probabilities of groups:
##    Lanjut  Drop Out    Pindah 
## 0.3647059 0.3294118 0.3058824 
## 
## Group means:
##               IPK     Skor
## Lanjut   3.403871 561.2258
## Drop Out 2.482500 447.0714
## Pindah   2.992692 446.2308
## 
## Coefficients of linear discriminants:
##               LD1         LD2
## IPK  -5.008766354  1.87668220
## Skor -0.008568593 -0.01445106
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9673 0.0327

# Koefisien asli dan terstandardisasi
koef_std <- lda_fit$scaling * sd_pooled(X, g)      
kable(data.frame(
  Peubah = rownames(lda_fit$scaling),
  `Asli LD1` = lda_fit$scaling[, 1], `Asli LD2` = lda_fit$scaling[, 2],
  `Baku LD1` = koef_std[, 1],        `Baku LD2` = koef_std[, 2],
  check.names = FALSE, row.names = NULL),
  digits = 4, caption = "Koefisien fungsi diskriminan LDA (asli dan terstandardisasi)")

Koefisien fungsi diskriminan LDA (asli dan terstandardisasi)
Peubah	Asli LD1	Asli LD2	Baku LD1	Baku LD2
IPK	-5.0088	1.8767	-0.9512	0.3564
Skor	-0.0086	-0.0145	-0.5181	-0.8738

# Persentase pemisahan yang dijelaskan tiap fungsi diskriminan
prop_trace <- lda_fit$svd^2 / sum(lda_fit$svd^2)

# Sentroid: rata-rata skor LD1 dan LD2 untuk tiap status
skor_lda <- data.frame(predict(lda_fit)$x, Status = g)
sentroid <- aggregate(cbind(LD1, LD2) ~ Status, data = skor_lda, FUN = mean)
kable(sentroid, digits = 3, caption = "Sentroid kelompok pada fungsi diskriminan")

Sentroid kelompok pada fungsi diskriminan
Status	LD1	LD2
Lanjut	-2.774	-0.246
Drop Out	2.819	-0.326
Pindah	0.271	0.644

Prior: peluang awal tiap status sama dengan proporsi sampel (0,365 untuk Lanjut, 0,329 Drop Out, 0,306 Pindah).
Kekuatan fungsi: LD1 menjelaskan 96,7% pemisahan antarkelompok, LD2 hanya 3,3%. Jadi hampir seluruh pembeda status tercakup dalam fungsi pertama.
Kontribusi peubah (koefisien baku): pada LD1, IPK (-0,951) dan skor tes (-0,518) sama-sama berkontribusi, dengan IPK lebih besar. Pada LD2, skor tes (-0,874) lebih berperan daripada IPK (0,356). Tanda negatif/positif hanya menunjukkan arah sumbu dan tidak mengandung arti baik atau buruk; yang dibaca adalah besar nilai mutlaknya.
Sentroid: nilai LD1 terpisah jelas antarstatus (lihat tabel sentroid), dengan Lanjut dan Drop Out berada di dua ujung yang berlawanan dan Pindah di antaranya.

ggplot(skor_lda, aes(x = LD1, y = LD2, colour = Status, shape = Status)) +
  geom_point(size = 2.6, alpha = 0.85) +
  stat_ellipse(level = 0.90, linewidth = 0.7) +
  geom_point(data = sentroid, aes(x = LD1, y = LD2), size = 5, shape = 4, stroke = 1.6,
             colour = "black", inherit.aes = FALSE) +               
  scale_colour_manual(values = warna) +
  labs(title = "Mahasiswa pada ruang diskriminan (LD1 vs LD2)",
       subtitle = "Tanda X hitam = sentroid kelompok") +
  theme_minimal(base_size = 12)

Sepanjang sumbu LD1 (horizontal), ketiga status tersusun berurutan dan Drop Out serta Lanjut berada di kedua ujung. Pindah berada di tengah dan tumpang tindih dengan Drop Out, sejalan dengan temuan Tukey bahwa skor tes tidak membedakan keduanya. Titik yang berada di area kelompok lain adalah calon kesalahan klasifikasi.

Model QDA

qda() dipanggil dengan cara yang sama seperti lda(). Bedanya, QDA menaksir matriks kovarians terpisah untuk tiap kelompok, sehingga batas keputusannya berupa kurva dan tidak ada koefisien tunggal yang bisa dibaca seperti pada LDA. Karena Box’s M menolak kesamaan ragam, QDA menjadi model utama.

qda_fit <- qda(Status ~ IPK + Skor, data = dat)   
qda_fit

## Call:
## qda(Status ~ IPK + Skor, data = dat)
## 
## Prior probabilities of groups:
##    Lanjut  Drop Out    Pindah 
## 0.3647059 0.3294118 0.3058824 
## 
## Group means:
##               IPK     Skor
## Lanjut   3.403871 561.2258
## Drop Out 2.482500 447.0714
## Pindah   2.992692 446.2308

Prior dan rata-rata kelompok sama dengan LDA. Perbedaannya ada pada cara model memperlakukan sebaran tiap kelompok: QDA mengizinkan setiap status memiliki bentuk dan ukuran sebaran sendiri. Karena tidak ada koefisien tunggal, kualitas QDA dinilai lewat hasil klasifikasi pada bagian berikut.

Evaluasi Kemampuan Klasifikasi

Kemampuan klasifikasi dinilai dengan dua cara.

Resubstitusi: model mengklasifikasikan ulang data yang dipakai membangunnya (predict(model)). Hasilnya cenderung terlalu optimistis.
LOOCV (leave-one-out cross-validation): setiap mahasiswa diklasifikasikan oleh model yang dibangun tanpa dirinya (CV = TRUE). Dengan 85 data yang terbatas, LOOCV lebih stabil daripada membagi data 70:30 dan lebih jujur daripada resubstitusi.

Ukuran yang dilaporkan: tingkat kesalahan = mean(prediksi != aktual), akurasi = mean(prediksi == aktual), dan sensitivitas tiap status = diag(tabel) / rowSums(tabel).

# Resubstitusi
lda_re <- predict(lda_fit)$class
qda_re <- predict(qda_fit)$class

# LOOCV (model dibangun ulang 85 kali, tiap kali tanpa satu mahasiswa)
lda_cv <- lda(Status ~ IPK + Skor, data = dat, CV = TRUE)$class
qda_cv <- qda(Status ~ IPK + Skor, data = dat, CV = TRUE)$class

ringkas_err <- data.frame(
  Model = c("LDA", "QDA"),
  `Salah (resubstitusi)` = c(sum(lda_re != g), sum(qda_re != g)),
  `Kesalahan (resubstitusi)` = c(mean(lda_re != g), mean(qda_re != g)),
  `Salah (LOOCV)` = c(sum(lda_cv != g), sum(qda_cv != g)),
  `Kesalahan (LOOCV)` = c(mean(lda_cv != g), mean(qda_cv != g)),
  `Akurasi (LOOCV)` = c(mean(lda_cv == g), mean(qda_cv == g)),
  check.names = FALSE)
kable_d(ringkas_err, c(3, 5, 6), 4, caption = "Kemampuan klasifikasi LDA dan QDA")

Kemampuan klasifikasi LDA dan QDA
Model	Salah (resubstitusi)	Kesalahan (resubstitusi)	Salah (LOOCV)	Kesalahan (LOOCV)	Akurasi (LOOCV)
LDA	7	0,0824	8	0,0941	0,9059
QDA	3	0,0353	4	0,0471	0,9529

kable(tabel_konfusi(g, lda_cv), caption = "Matriks klasifikasi LDA (LOOCV)")

Matriks klasifikasi LDA (LOOCV)
Aktual ↓ / Prediksi →	Lanjut	Drop Out	Pindah
Lanjut	27	0	4
Drop Out	0	26	2
Pindah	1	1	24

kable(tabel_konfusi(g, qda_cv), caption = "Matriks klasifikasi QDA (LOOCV)")

Matriks klasifikasi QDA (LOOCV)
Aktual ↓ / Prediksi →	Lanjut	Drop Out	Pindah
Lanjut	30	0	1
Drop Out	0	27	1
Pindah	1	1	24

# Sensitivitas: persentase tiap status yang dikenali dengan benar (LOOCV)
sens <- data.frame(
  Status = levels(g),
  `Sensitivitas LDA` = as.vector(diag(table(g, lda_cv)) / table(g)),
  `Sensitivitas QDA` = as.vector(diag(table(g, qda_cv)) / table(g)),
  check.names = FALSE)
kable_d(sens, 2:3, 3, caption = "Sensitivitas per status (LOOCV)")

Sensitivitas per status (LOOCV)
Status	Sensitivitas LDA	Sensitivitas QDA
Lanjut	0,871	0,968
Drop Out	0,929	0,964
Pindah	0,923	0,923

LOOCV: kesalahan LDA = 0,0941 (8 dari 85 salah) dan QDA = 0,0471 (4 dari 85 salah). QDA lebih rendah kesalahannya, konsisten dengan hasil Box’s M bahwa model kuadratik lebih sesuai untuk data ini.
Resubstitusi: kesalahan LDA 0,0824 dan QDA 0,0353. Angka ini lebih kecil daripada LOOCV (sifat optimistis resubstitusi), sehingga tidak dijadikan dasar utama.
Pola kesalahan: pada matriks klasifikasi, kesalahan LDA terbanyak adalah mahasiswa Lanjut yang diprediksi Pindah (4 kasus), hal ini terjadi karena LDA memakai satu matriks ragam gabungan yang lebih sempit daripada ragam skor kelompok Lanjut (lihat tabel komponen pada Box’s M), sehingga sebagian mahasiswa Lanjut yang skornya menyebar jauh dari pusat kelompok “terlempar” ke kelompok lain. Pada QDA kesalahan menyebar tipis, hanya satu atau dua kasus per sel.
Sensitivitas: status dengan sensitivitas QDA terendah adalah Pindah (92,3%). Ini wajar karena kelompok Pindah tumpang tindih dengan Drop Out pada skor tes.

Apakah klasifikasi lebih baik daripada menebak?

Dua kriteria baku menilai apakah akurasi model benar-benar bermakna dibanding tebakan acak:

Kriteria peluang proporsional (C_pro = Σ pₖ²): akurasi yang diperoleh dengan menebak sesuai proporsi kelompok. Akurasi model dianggap baik bila minimal 1,25 × C_pro.
Statistik Q Press = (N − n·K)² / (N(K − 1)), dengan n jumlah klasifikasi benar dan K jumlah kelompok. Dibandingkan dengan khi-kuadrat db = 1 (nilai kritis 3,84 pada taraf 5%; 6,63 pada 1%).

N        <- nrow(dat)
c_pro    <- sum((n_kelompok / N)^2)                                  
benar    <- c(LDA = sum(lda_cv == g), QDA = sum(qda_cv == g))         
press_q  <- (N - benar * k)^2 / (N * (k - 1))                        

kable(data.frame(
  Model = names(benar),
  `Akurasi LOOCV` = benar / N,
  `Batas 1,25 x C_pro` = 1.25 * c_pro,
  `Q Press` = press_q,
  `p-value` = fmt_p(pchisq(press_q, df = 1, lower.tail = FALSE)),
  check.names = FALSE, row.names = NULL),
  digits = 3, caption = "Perbandingan akurasi dengan tebakan acak")

Perbandingan akurasi dengan tebakan acak
Model	Akurasi LOOCV	Batas 1,25 x C_pro	Q Press	p-value
LDA	0.906	0.419	125.388	< 0,0001
QDA	0.953	0.419	146.847	< 0,0001

Menebak sesuai proporsi kelompok hanya menghasilkan akurasi sekitar 33,5% (C_pro = 0,335), dan batas 1,25 × C_pro adalah 41,9%. Akurasi LOOCV model (90,6% untuk LDA dan 95,3% untuk QDA) jauh di atas batas tersebut, dan Q Press (125,4 dan 146,8) jauh melampaui nilai kritis 6,63. Jadi kemampuan klasifikasi kedua model bukan kebetulan.

Batas Keputusan

expand.grid() membuat ribuan kombinasi IPK dan skor tes yang menutupi seluruh bidang. Model memprediksi status untuk tiap kombinasi (predict(model, grid)$class), lalu geom_raster() mewarnai bidang menurut prediksi. Titik amatan ditumpangkan di atas warna. Dua panel (facet_wrap) membandingkan batas lurus LDA dengan batas melengkung QDA.

grid <- expand.grid(
  IPK  = seq(min(dat$IPK) - 0.1, max(dat$IPK) + 0.1, length.out = 200),
  Skor = seq(min(dat$Skor) - 20, max(dat$Skor) + 20, length.out = 200)
)

grid_semua <- rbind(
  transform(grid, Prediksi = predict(lda_fit, grid)$class, Model = "LDA"),   
  transform(grid, Prediksi = predict(qda_fit, grid)$class, Model = "QDA")    
)

ggplot() +
  geom_raster(data = grid_semua, aes(x = IPK, y = Skor, fill = Prediksi), alpha = 0.28) +
  geom_point(data = dat, aes(x = IPK, y = Skor, colour = Status, shape = Status), size = 2.2) +
  scale_fill_manual(values = warna, name = "Wilayah klasifikasi") +
  scale_colour_manual(values = warna) +
  facet_wrap(~ Model) +
  labs(title = "Batas keputusan pada bidang IPK dan skor tes",
       x = "IPK S1", y = "Skor tes masuk S2") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")

Wilayah Lanjut (hijau) menempati sisi kanan-atas: semakin tinggi IPK, semakin rendah skor tes yang dibutuhkan, dan sebaliknya. Artinya kedua peubah saling mengompensasi. Pada LDA, batas antarwilayah berupa garis lurus; pada QDA, batasnya melengkung mengikuti bentuk sebaran tiap kelompok, terutama di sekitar wilayah Pindah dan Drop Out. Titik-titik yang berada di wilayah warna lain adalah mahasiswa yang salah diklasifikasikan.

Penyusunan Kriteria Seleksi

Bagian ini menerjemahkan model menjadi aturan yang dapat dipakai Ketua Program Studi. Fokusnya adalah mengenali calon yang berpeluang Lanjut, karena Drop Out dan Pindah tidak dapat dibedakan dengan baik oleh kedua peubah ini.

Peluang Lanjut untuk calon mahasiswa baru

predict(qda_fit, newdata = kandidat) memberikan dua hal untuk setiap calon: kelas prediksi ($class) dan peluang posterior tiap status ($posterior). Peluang posterior menyatakan seberapa yakin model bahwa calon tersebut akan berstatus Lanjut, Drop Out, atau Pindah. Enam calon fiktif dengan kombinasi IPK dan skor berbeda dipakai sebagai contoh.

kandidat <- data.frame(
  Calon = LETTERS[1:6],
  IPK   = c(3.6, 3.3, 3.0, 2.8, 3.4, 2.5),
  Skor  = c(450, 500, 560, 600, 400, 500)
)

pred_kand <- predict(qda_fit, newdata = kandidat)     
hasil_kand <- data.frame(
  kandidat,
  `P(Lanjut)`   = pred_kand$posterior[, "Lanjut"],
  `P(Drop Out)` = pred_kand$posterior[, "Drop Out"],
  `P(Pindah)`   = pred_kand$posterior[, "Pindah"],
  Prediksi      = as.character(pred_kand$class),
  check.names = FALSE)
tampil_kand <- hasil_kand                                  
tampil_kand$IPK <- id(tampil_kand$IPK, 1)                  
tampil_kand[4:6] <- lapply(tampil_kand[4:6], id, d = 3)   
kable(tampil_kand, caption = "Peluang posterior QDA untuk enam calon contoh")

Peluang posterior QDA untuk enam calon contoh
Calon	IPK	Skor	P(Lanjut)	P(Drop Out)	P(Pindah)	Prediksi
A	3,6	450	1,000	0,000	0,000	Lanjut
B	3,3	500	0,994	0,000	0,006	Lanjut
C	3,0	560	0,925	0,014	0,062	Lanjut
D	2,8	600	0,527	0,312	0,162	Lanjut
E	3,4	400	0,407	0,000	0,593	Pindah
F	2,5	500	0,000	0,975	0,025	Drop Out

Calon yang diprediksi Lanjut: A, B, C, D. Beberapa hal penting dari tabel:

Calon A (IPK 3,6, skor 450) hampir pasti Lanjut walaupun skornya tidak tinggi, karena IPK-nya tinggi. Sebaliknya calon E (IPK 3,4, skor 400) diprediksi Pindah karena skornya rendah, walaupun IPK-nya baik. Ini menunjukkan kedua peubah harus dilihat bersamaan.
Calon D (IPK 2,8, skor 600) memiliki peluang Lanjut hanya 0,527, yaitu mendekati batas. IPK 2,8 juga di bawah IPK terendah mahasiswa Lanjut pada data (2,96), sehingga prediksinya berada di luar jangkauan data dan sebaiknya tidak dipercaya penuh.
Calon F (IPK 2,5) diprediksi Drop Out dengan peluang tinggi.

Ambang skor tes minimum menurut IPK

Untuk tiap nilai IPK dari 3,0 sampai 3,8, dicari skor tes terendah yang membuat peluang Lanjut mencapai ambang tertentu (0,5; 0,7; dan 0,9).

expand.grid(IPK = ..., Skor = 250:800) membuat kombinasi IPK dan skor satu per satu angka.
predict(qda_fit, grid)$posterior[, "Lanjut"] menghitung peluang Lanjut tiap kombinasi.
Fungsi min_skor() mengambil skor terkecil yang peluangnya mencapai ambang.
IPK di bawah 3,0 tidak dimasukkan karena hampir tidak ada mahasiswa Lanjut di wilayah itu, sehingga model akan mengekstrapolasi.

ipk_seq <- seq(3.0, 3.8, by = 0.1)                                  
g_ambang <- expand.grid(IPK = ipk_seq, Skor = 250:800)               
g_ambang$P <- predict(qda_fit, g_ambang)$posterior[, "Lanjut"]       

min_skor <- function(ipk, ambang) {
  d <- g_ambang[abs(g_ambang$IPK - ipk) < 1e-9 & g_ambang$P >= ambang, ]    
  if (nrow(d) == 0) NA else min(d$Skor)                                   
}

tabel_ambang <- data.frame(
  IPK = ipk_seq,
  `P(Lanjut) >= 0,5` = sapply(ipk_seq, min_skor, ambang = 0.5),
  `P(Lanjut) >= 0,7` = sapply(ipk_seq, min_skor, ambang = 0.7),
  `P(Lanjut) >= 0,9` = sapply(ipk_seq, min_skor, ambang = 0.9),
  check.names = FALSE)
kable(tabel_ambang, digits = c(1, 0, 0, 0),
      caption = "Skor tes minimum agar peluang Lanjut mencapai ambang, menurut IPK")

Skor tes minimum agar peluang Lanjut mencapai ambang, menurut IPK
IPK	P(Lanjut) >= 0,5	P(Lanjut) >= 0,7	P(Lanjut) >= 0,9
3.0	528	540	557
3.1	498	510	526
3.2	468	480	497
3.3	438	450	468
3.4	407	420	438
3.5	376	390	409
3.6	344	359	379
3.7	312	328	349
3.8	279	296	318

Tabel menunjukkan tukar-tambah (trade-off) antara IPK dan skor tes. Misalnya, agar peluang Lanjut minimal 0,5, mahasiswa ber-IPK 3,0 membutuhkan skor tes sekitar 528, sedangkan yang ber-IPK 3,5 cukup dengan sekitar 376. Setiap kenaikan IPK 0,1 menurunkan skor tes yang diperlukan sekitar 31 poin. Semakin ketat ambang (0,9), semakin tinggi skor yang diminta pada IPK yang sama.

g_kontur <- expand.grid(IPK = seq(2.0, 4.0, length.out = 200), Skor = seq(250, 750, length.out = 200))
g_kontur$P <- predict(qda_fit, g_kontur)$posterior[, "Lanjut"]       

ggplot() +
  geom_point(data = dat, aes(x = IPK, y = Skor, colour = Status, shape = Status), size = 2.4, alpha = 0.85) +
  geom_contour(data = g_kontur, aes(x = IPK, y = Skor, z = P), breaks = 0.5, colour = "black", linewidth = 0.9) +
  geom_contour(data = g_kontur, aes(x = IPK, y = Skor, z = P), breaks = 0.7, colour = "black", linewidth = 0.7, linetype = "dashed") +
  geom_contour(data = g_kontur, aes(x = IPK, y = Skor, z = P), breaks = 0.9, colour = "black", linewidth = 0.7, linetype = "dotted") +
  scale_colour_manual(values = warna) +
  labs(title = "Kriteria seleksi berdasarkan peluang Lanjut (QDA)",
       subtitle = "Calon di sisi kanan-atas kurva memiliki peluang Lanjut lebih tinggi",
       x = "IPK S1", y = "Skor tes masuk S2") +
  theme_minimal(base_size = 12)

Pada rentang IPK 3,0 sampai 3,8 (rentang data mahasiswa Lanjut), kurva tegas (peluang 0,5) hampir lurus dan memisahkan titik hijau (Lanjut) dari titik oranye dan ungu. Hampir semua mahasiswa Lanjut berada di sisi kanan-atas kurva, sedangkan hampir semua Drop Out dan Pindah berada di sisi kiri-bawah. Kurva 0,7 dan 0,9 bergeser ke kanan-atas: makin ketat kriteria, makin sedikit calon yang lolos tetapi makin yakin bahwa mereka akan Lanjut.

Aturan praktis yang mudah dipakai

Kurva ambang 0,5 hampir lurus, sehingga dapat dirangkum dalam satu rumus sederhana. lm() (regresi linear) pada tabel ambang menghasilkan intersep a dan kemiringan b, sehingga aturannya menjadi Skor ≥ a + b × IPK. Aturan ini kemudian diuji pada 85 mahasiswa untuk melihat seberapa baik ia memisahkan Lanjut dari bukan Lanjut (Drop Out atau Pindah). ifelse() memberi label prediksi dan table() menghitung ketepatannya.

# Aturan: calon "Diterima" bila Skor >= a + b x IPK
dat$Aturan <- ifelse(dat$Skor >= a_amb + b_amb * dat$IPK, "Diterima", "Tidak diterima")
dat$Hasil  <- ifelse(dat$Status == "Lanjut", "Lanjut", "Bukan Lanjut")

tab_aturan <- table(`Status sebenarnya` = dat$Hasil, `Keputusan aturan` = dat$Aturan)
kable(as.data.frame.matrix(tab_aturan), caption = "Ketepatan aturan praktis terhadap hasil sebenarnya")

Ketepatan aturan praktis terhadap hasil sebenarnya
	Diterima	Tidak diterima
Bukan Lanjut	1	53
Lanjut	30	1

# Ukuran ketepatan aturan
sens_a <- tab_aturan["Lanjut", "Diterima"] / sum(tab_aturan["Lanjut", ])            
spes_a <- tab_aturan["Bukan Lanjut", "Tidak diterima"] / sum(tab_aturan["Bukan Lanjut", ])  
akur_a <- sum(diag(tab_aturan[c("Lanjut", "Bukan Lanjut"), c("Diterima", "Tidak diterima")])) / N

Aturan praktis hasil pendekatan adalah:

Calon diterima bila Skor tes ≥ 1462 − 311 × IPK

Pada 85 mahasiswa, aturan ini menerima 96,8% mahasiswa yang memang Lanjut (sensitivitas), menolak 98,1% mahasiswa yang bukan Lanjut (spesifisitas), dengan akurasi keseluruhan 97,6%. Karena aturan ini dirumuskan dari data yang sama, angkanya bersifat optimistis (tidak divalidasi silang), tetapi ia cukup mewakili batas yang dibentuk QDA dan jauh lebih mudah dipakai dalam seleksi.

Ambang 0,5 menyeimbangkan dua jenis kesalahan: menolak calon yang sebenarnya akan Lanjut dan menerima calon yang sebenarnya tidak Lanjut. Bila program studi ingin sangat selektif (lebih mementingkan tidak menerima calon berisiko), pakai ambang 0,7 atau 0,9 pada tabel ambang. Bila daya tampung besar dan ingin lebih longgar, pakai ambang 0,5. Pilihan ambang adalah keputusan kebijakan, bukan keputusan statistik.

Kesimpulan

kable(data.frame(
  Aspek = c("Normalitas multivariat (Mardia)", "Kesamaan ragam-peragam (Box's M)",
            "Perbedaan antarstatus (MANOVA)", "Kesalahan LOOCV LDA", "Kesalahan LOOCV QDA"),
  Hasil = c(ifelse(n_tolak_mardia == 0, "Terpenuhi", paste0("Tidak terpenuhi pada ", paste(grup_tolak, collapse = ", "))),
            paste0("Ditolak (p ", fmt_p(bm$p), ")"),
            paste0("Berbeda nyata (Wilks ", id(st[1, "Wilks"], 3), ", p ", fmt_p(st[1, "Pr(>F)"]), ")"),
            id(ringkas_err[1, 5], 4), id(ringkas_err[2, 5], 4))
), caption = "Ringkasan hasil analisis")

Ringkasan hasil analisis
Aspek	Hasil
Normalitas multivariat (Mardia)	Tidak terpenuhi pada Pindah
Kesamaan ragam-peragam (Box’s M)	Ditolak (p 0,0134)
Perbedaan antarstatus (MANOVA)	Berbeda nyata (Wilks 0,126, p < 0,0001)
Kesalahan LOOCV LDA	0,0941
Kesalahan LOOCV QDA	0,0471

Asumsi. Matriks ragam-peragam ketiga status tidak sama (Box’s M, p 0,0134), sehingga QDA lebih sesuai secara teori daripada LDA. Normalitas multivariat tidak terpenuhi pada kelompok Pindah (pelanggaran ringan, terkait beberapa pencilan).
Perbedaan status. Ketiga status berbeda nyata pada IPK dan skor tes secara bersama-sama (MANOVA, Wilks = 0,126), sehingga analisis diskriminan layak.
Peran peubah. IPK S1 membedakan ketiga status secara bertingkat dan menjadi pembeda keseluruhan yang lebih kuat. Skor tes masuk S2 terutama membedakan Lanjut dari Drop Out dan Pindah, tetapi tidak membedakan Drop Out dari Pindah. Fungsi diskriminan pertama menjelaskan 96,7% pemisahan.
Kinerja model. Tingkat kesalahan LOOCV adalah 0,0471 untuk QDA (akurasi 95,3%) dan 0,0941 untuk LDA (akurasi 90,6%), jauh di atas akurasi tebakan acak (33,5%). QDA lebih baik, sejalan dengan hasil Box’s M.
Kriteria seleksi. Gunakan kedua peubah bersamaan, karena keduanya saling mengompensasi. Aturan praktis: terima calon bila Skor tes ≥ 1462 − 311 × IPK (setara peluang Lanjut sekitar 0,5 menurut QDA). Untuk seleksi yang lebih ketat gunakan ambang 0,7 atau 0,9 pada tabel ambang. Kriteria sebaiknya difokuskan mengenali calon yang berpeluang Lanjut, bukan membedakan Pindah dari Drop Out, karena keduanya tidak dapat dipisahkan dengan baik oleh IPK dan skor tes.

Keterbatasan

Ukuran sampel kecil (85 mahasiswa; 26 sampai 31 per kelompok), sehingga taksiran matriks ragam-peragam pada QDA kurang stabil.
Normalitas multivariat tidak terpenuhi sepenuhnya pada kelompok Pindah, dan Box’s M menolak kesamaan ragam.
Kelompok Pindah dan Drop Out saling tumpang tindih, sehingga kemampuan klasifikasi keduanya terbatas.
Model hanya memakai dua peubah. Faktor lain seperti latar belakang program studi S1, motivasi, kebiasaan belajar, atau dukungan pembiayaan belum dimasukkan, padahal dapat memengaruhi status Pindah dan Drop Out.
Aturan praktis dan ambang peluang diturunkan dari data yang sama sehingga belum tervalidasi pada angkatan lain. Sebaiknya diuji ulang pada angkatan berikutnya.
Prediksi untuk IPK di luar jangkauan data (misalnya di bawah 3,0 pada kelompok Lanjut) bersifat ekstrapolasi dan tidak boleh dipercaya penuh.
Analisis ini bersifat prediktif/asosiatif, bukan sebab-akibat: IPK dan skor tes membantu memprediksi status, tetapi tidak membuktikan penyebab Drop Out atau Pindah.

Analisis Diskriminan Keberhasilan Studi Mahasiswa S2 Semester Pertama

Tugas Praktikum Teknik Peubah Ganda

Fadly Syabani

09-11-2026