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(tidyr)    # untuk mengubah bentuk data: pivot_longer()

## Warning: package 'tidyr' was built under R version 4.5.2

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 lebih singkat dan konsisten.

id() membulatkan angka lalu menampilkannya dengan koma desimal, dipakai pada teks interpretasi.
fmt_p() menampilkan p-value; nilai yang sangat kecil ditulis < 0,0001.
sd_pooled() menghitung simpangan baku gabungan dalam kelompok, dipakai 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 spesies 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 = kelas aktual, kolom = kelas 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 spesies
warna <- c(setosa = "#1b9e77", versicolor = "#d95f02", virginica = "#7570b3")

Impor Data

Memuat dataset iris

iris adalah dataset bawaan R, sehingga cukup dipanggil dengan data(iris). Salinannya disimpan sebagai iris_data agar data asli tidak berubah. Fungsi str() menampilkan struktur data (jumlah baris, nama kolom, tipe data), head() menampilkan enam baris pertama, dan is.na() dipakai untuk memeriksa data yang hilang.

data(iris)               
iris_data <- iris        

str(iris_data)

## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

kable(head(iris_data), caption = "Enam baris pertama dataset iris")

Enam baris pertama dataset iris
Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
5.1	3.5	1.4	0.2	setosa
4.9	3.0	1.4	0.2	setosa
4.7	3.2	1.3	0.2	setosa
4.6	3.1	1.5	0.2	setosa
5.0	3.6	1.4	0.2	setosa
5.4	3.9	1.7	0.4	setosa

# Jumlah data hilang (NA) di seluruh dataset
cat("Jumlah data hilang (NA):", sum(is.na(iris_data)), "\n")

## Jumlah data hilang (NA): 0

Data berisi 150 amatan dan 5 kolom. Empat kolom pertama bertipe numerik (num) dan menjadi peubah pengukuran (panjang dan lebar sepal serta petal, dalam cm). Kolom Species bertipe faktor dengan 3 level dan menjadi peubah kelompok yang akan diprediksi. Jumlah data hilang adalah 0, sehingga data langsung dapat dipakai tanpa pembersihan.

Ringkasan statistik

summary() menghitung nilai minimum, kuartil, median, rata-rata, dan maksimum tiap peubah. table() menghitung banyaknya amatan per spesies untuk memeriksa keseimbangan kelompok.

kable(summary(iris_data[, 1:4]), caption = "Ringkasan statistik deskriptif empat peubah")

Ringkasan statistik deskriptif empat peubah
Sepal.Length	Sepal.Width	Petal.Length	Petal.Width
Min. :4.300	Min. :2.000	Min. :1.000	Min. :0.100
1st Qu.:5.100	1st Qu.:2.800	1st Qu.:1.600	1st Qu.:0.300
Median :5.800	Median :3.000	Median :4.350	Median :1.300
Mean :5.843	Mean :3.057	Mean :3.758	Mean :1.199
3rd Qu.:6.400	3rd Qu.:3.300	3rd Qu.:5.100	3rd Qu.:1.800
Max. :7.900	Max. :4.400	Max. :6.900	Max. :2.500

kable(as.data.frame(table(iris_data$Species)),
      col.names = c("Spesies", "Jumlah"),
      caption = "Jumlah amatan per spesies")

Jumlah amatan per spesies
Spesies	Jumlah
setosa	50
versicolor	50
virginica	50

Rata-rata panjang sepal adalah 5,84 cm, sedangkan rata-rata panjang petal 3,76 cm dengan rentang 1,0 sampai 6,9 cm. Rentang petal yang lebar menandakan ukuran petal sangat bervariasi antarbunga, kemungkinan karena perbedaan spesies. Setiap spesies memiliki 50 amatan, sehingga data seimbang dan tidak ada spesies yang mendominasi model.

Eksplorasi Data

Rata-rata dan simpangan baku per spesies

group_by(Species) memisahkan data menurut spesies, lalu summarise(across(...)) menghitung statistik untuk seluruh kolom numerik sekaligus. Hasilnya dibulatkan 3 desimal dengan round().

rata2 <- iris_data %>%
  group_by(Species) %>%                                          
  summarise(across(everything(), ~ round(mean(.x), 3)), .groups = "drop")  

simpangan <- iris_data %>%
  group_by(Species) %>%
  summarise(across(everything(), ~ round(sd(.x), 3)), .groups = "drop")    

kable(rata2, caption = "Rata-rata peubah per spesies")

Rata-rata peubah per spesies
Species	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width
setosa	5.006	3.428	1.462	0.246
versicolor	5.936	2.770	4.260	1.326
virginica	6.588	2.974	5.552	2.026

kable(simpangan, caption = "Simpangan baku peubah per spesies")

Simpangan baku peubah per spesies
Species	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width
setosa	0.352	0.379	0.174	0.105
versicolor	0.516	0.314	0.470	0.198
virginica	0.636	0.322	0.552	0.275

Setosa memiliki petal paling kecil (panjang 1,46 cm, lebar 0,25 cm), tetapi sepal paling lebar (3,43 cm).
Virginica memiliki petal paling besar (panjang 5,55 cm, lebar 2,03 cm).
Versicolor berada di antara keduanya pada hampir semua peubah.
Simpangan baku petal setosa jauh lebih kecil dibanding spesies lain. Perbedaan keragaman antarspesies ini akan muncul kembali pada uji kesamaan matriks ragam-peragam.

Sebaran peubah menurut spesies

pivot_longer(-Species) mengubah data dari bentuk lebar (4 kolom peubah) menjadi bentuk panjang (satu kolom peubah, satu kolom nilai), yang dibutuhkan ggplot2 untuk membuat panel per peubah dengan facet_wrap(). geom_boxplot() menggambar kotak sebaran dan geom_jitter() menambahkan titik amatan sebenarnya.

iris_data %>%
  pivot_longer(-Species, names_to = "peubah", values_to = "nilai") %>%   
  ggplot(aes(x = Species, y = nilai, fill = Species)) +
  geom_boxplot(alpha = 0.65, outlier.shape = NA, show.legend = FALSE) + 
  geom_jitter(width = 0.12, height = 0, size = 1, alpha = 0.5, show.legend = FALSE) + 
  scale_fill_manual(values = warna) +
  facet_wrap(~ peubah, scales = "free_y") +                              
  labs(title = "Sebaran peubah iris menurut spesies", x = NULL, y = "Nilai (cm)") +
  theme_minimal(base_size = 11)

Pada panjang dan lebar petal, kotak setosa terpisah sepenuhnya dari dua spesies lain, sehingga petal sangat efektif untuk mengenali setosa. Versicolor dan virginica juga cukup terpisah pada petal, tetapi masih ada sedikit daerah tumpang tindih. Pada lebar sepal, tumpang tindih versicolor dan virginica lebih besar, sehingga peubah ini kurang membedakan keduanya. Secara visual, peubah petal memberi pemisahan terkuat.

Hubungan antarpeubah

Data sepal dan petal digabung menjadi satu data berlabel panel agar kedua diagram pencar tampil berdampingan. stat_ellipse() menggambar elips yang melingkupi sekitar 95% sebaran tiap spesies. cor() menghitung korelasi Pearson antarpeubah.

sepal <- data.frame(x = iris_data$Sepal.Length, y = iris_data$Sepal.Width,
                    Species = iris_data$Species, panel = "Sepal (panjang vs lebar)")
petal <- data.frame(x = iris_data$Petal.Length, y = iris_data$Petal.Width,
                    Species = iris_data$Species, panel = "Petal (panjang vs lebar)")

ggplot(rbind(sepal, petal), aes(x = x, y = y, colour = Species)) +
  geom_point(size = 2, alpha = 0.8) +
  stat_ellipse(level = 0.95, linewidth = 0.7) +      
  scale_colour_manual(values = warna) +
  facet_wrap(~ panel, scales = "free") +
  labs(title = "Diagram pencar menurut spesies", x = "Panjang (cm)", y = "Lebar (cm)") +
  theme_minimal(base_size = 11)

kable(round(cor(iris_data[, 1:4]), 3), caption = "Matriks korelasi antarpeubah (seluruh data)")

Matriks korelasi antarpeubah (seluruh data)
	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width
Sepal.Length	1.000	-0.118	0.872	0.818
Sepal.Width	-0.118	1.000	-0.428	-0.366
Petal.Length	0.872	-0.428	1.000	0.963
Petal.Width	0.818	-0.366	0.963	1.000

Pada diagram petal, ketiga spesies membentuk tiga gerombol yang hampir terpisah sempurna; hanya versicolor dan virginica yang sedikit bersinggungan. Pada diagram sepal, gerombol saling tumpang tindih lebih banyak. Ukuran dan bentuk elips juga berbeda antarspesies, yang kembali mengisyaratkan matriks ragam-peragam yang tidak sama. Korelasi panjang dan lebar petal sangat tinggi (r = 0,963). Pada analisis diskriminan hal ini tidak merusak ketepatan klasifikasi, tetapi koefisien kedua peubah sulit ditafsirkan secara terpisah.

Pembagian Data Latih dan Data Uji

Data dibagi 70% latih dan 30% uji secara berstrata, yaitu pembagian dilakukan di dalam tiap spesies agar proporsi spesies sama di kedua himpunan.

set.seed(2026) mengunci bilangan acak sehingga hasil dapat direproduksi.
split(seq_len(nrow(iris_data)), iris_data$Species) mengelompokkan nomor baris menurut spesies.
sample(ii, size = round(0.7 * length(ii))) mengambil 70% nomor baris secara acak di tiap spesies.
unlist() menggabungkan hasilnya menjadi satu vektor nomor baris data latih.
Baris yang tidak terpilih (-idx_latih) menjadi data uji.

set.seed(2026)                                 

idx_latih <- unlist(lapply(
  split(seq_len(nrow(iris_data)), iris_data$Species),         
  function(ii) sample(ii, size = round(0.7 * length(ii)))     
))
idx_latih <- sort(idx_latih)                     

latih <- iris_data[idx_latih, ]                 
uji   <- iris_data[-idx_latih, ]                 

kable(data.frame(
  Spesies = levels(iris_data$Species),
  Latih   = as.vector(table(latih$Species)),
  Uji     = as.vector(table(uji$Species))
), caption = "Jumlah amatan data latih dan data uji per spesies")

Jumlah amatan data latih dan data uji per spesies
Spesies	Latih	Uji
setosa	35	15
versicolor	35	15
virginica	35	15

Data latih berisi 105 amatan (35 per spesies) dan data uji berisi 45 amatan (15 per spesies). Karena pembagian berstrata, proporsi spesies sama di kedua himpunan, sehingga perbandingan kinerja model tidak bias oleh perbedaan jumlah spesies. Data uji tidak boleh dipakai untuk membangun model maupun menguji asumsi sehingga disimpan khusus untuk evaluasi akhir.

Uji Asumsi pada Data Latih

Dua asumsi utama analisis diskriminan diperiksa hanya pada data latih. Matriks X_latih menyimpan keempat peubah, sedangkan g_latih menyimpan label spesies.

X_latih <- as.matrix(latih[, 1:4])   
g_latih <- latih$Species

Uji normalitas multivariat (Mardia)

Uji Mardia memeriksa dua hal pada tiap spesies:

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 spesies dengan lapply() dan hasilnya digabung dengan rbind().

Hipotesis: H₀ = data spesies 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 spesies lalu gabungkan
hasil_mardia <- do.call(rbind, lapply(levels(g_latih), function(l) {
  cbind(Spesies = l, mardia(X_latih[g_latih == 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("Spesies", "n", "Skewness (b1)", "p Skewness", "Kurtosis (b2)",
                    "z Kurtosis", "p Kurtosis", "Keputusan"),
      caption = "Uji normalitas multivariat Mardia per spesies (data latih)")

Uji normalitas multivariat Mardia per spesies (data latih)
Spesies	n	Skewness (b1)	p Skewness	Kurtosis (b2)	z Kurtosis	p Kurtosis	Keputusan
setosa	35	3.625	0.3885	25.128	0.481	0.6302	Gagal tolak H0
versicolor	35	3.761	0.3440	22.465	-0.655	0.5122	Gagal tolak H0
virginica	35	3.589	0.4007	23.479	-0.223	0.8238	Gagal tolak H0

Nilai p terkecil di seluruh uji adalah 0,3440. Seluruh p-value berada di atas 0,05, sehingga H₀ tidak ditolak pada ketiga spesies: asumsi normal ganda dianggap terpenuhi. Dengan hanya 35 amatan per spesies, daya uji Mardia terbatas, jadi hasil “tidak ditolak” berarti tidak ada bukti kuat pelanggaran, bukan bukti bahwa data benar-benar normal.

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

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

Si adalah daftar matriks kovarians tiap spesies, Sp adalah matriks kovarians gabungan (pooled).
M membandingkan determinan matriks gabungan dengan determinan tiap matriks spesies. 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 spesies 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),
       logdet = setNames(sapply(Si, function(S) log(det(S))), gs))
}

bm <- box_m(X_latih, g_latih)

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 (data latih)")

Uji Box’s M untuk kesamaan matriks ragam-peragam (data latih)
Statistik	Nilai
M	117,225
Khi-kuadrat terkoreksi	110,636
Derajat bebas	20
p-value	< 0,0001

# Log-determinan matriks kovarians tiap spesies (ukuran 'volume' sebaran)
kable(data.frame(Spesies = names(bm$logdet), `Log determinan` = round(bm$logdet, 3),
                 check.names = FALSE, row.names = NULL),
      caption = "Log-determinan matriks kovarians per spesies")

Log-determinan matriks kovarians per spesies
Spesies	Log determinan
setosa	-13.069
versicolor	-10.798
virginica	-8.751

Statistik Box’s M terkoreksi bernilai 110,636 dengan derajat bebas 20 dan p-value < 0,0001. Karena p-value < 0,05, H₀ ditolak: matriks ragam-peragam antarspesies tidak sama. Tabel log-determinan memperlihatkan perbedaannya, yaitu volume sebaran tiap spesies tidak seragam. Akibatnya asumsi dasar LDA tidak terpenuhi dan QDA lebih sesuai secara teori; LDA tetap dibangun sebagai pembanding, dan pilihan akhir akan dinilai dari kinerja prediksi.

Uji Box’s M sangat peka terhadap pelanggaran normalitas dan cenderung mudah menolak H₀. Karena itu keputusan memilih LDA atau QDA tidak berhenti pada uji ini saja; kinerja prediksi pada data uji tetap menjadi penentu

Uji perbedaan vektor rata-rata (MANOVA)

Analisis diskriminan hanya bermakna jika rata-rata kelompok memang berbeda. manova() menguji H₀: vektor rata-rata keempat peubah sama di ketiga spesies. summary(..., test = "Wilks") memakai statistik Wilks’ Lambda; nilainya mendekati 0 bila kelompok sangat terpisah dan mendekati 1 bila tidak ada pemisahan.

fit_manova <- manova(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species,
                     data = latih)
hasil_manova <- summary(fit_manova, test = "Wilks")
st <- hasil_manova$stats                      

kable(data.frame(
  Sumber        = c("Species", "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) pada data latih")

MANOVA (uji Wilks’ Lambda) pada data latih
Sumber	Derajat bebas	Wilks’ Lambda	F hampiran	df1; df2	p-value
Species	2	0,0256	130,059	8; 198	< 0,0001
Galat (Residuals)	102

Wilks’ Lambda bernilai 0,0256 dengan p-value < 0,0001. Nilai yang sangat dekat dengan 0 dan p yang jauh di bawah 0,05 berarti vektor rata-rata ketiga spesies berbeda nyata, sehingga ada dasar yang kuat untuk membangun fungsi diskriminan.

Ringkasan uji asumsi

kable(data.frame(
  Asumsi = c("Normal ganda (Mardia)", "Kesamaan matriks ragam-peragam (Box's M)",
             "Perbedaan vektor rata-rata (MANOVA)"),
  Hasil = c(paste0("Tolak H0 pada ", n_tolak_mardia, " dari 3 spesies"),
            paste0("p ", fmt_p(bm$p)),
            paste0("p ", fmt_p(hasil_manova$stats[1, "Pr(>F)"]))),
  Kesimpulan = c(ifelse(n_tolak_mardia == 0, "Terpenuhi", "Tidak sepenuhnya terpenuhi"),
                 ifelse(bm$p < 0.05, "Tidak terpenuhi (pakai QDA)", "Terpenuhi (pakai LDA)"),
                 ifelse(hasil_manova$stats[1, "Pr(>F)"] < 0.05, "Kelompok berbeda nyata", "Tidak berbeda"))
), caption = "Ringkasan pemeriksaan asumsi")

Ringkasan pemeriksaan asumsi
Asumsi	Hasil	Kesimpulan
Normal ganda (Mardia)	Tolak H0 pada 0 dari 3 spesies	Terpenuhi
Kesamaan matriks ragam-peragam (Box’s M)	p < 0,0001	Tidak terpenuhi (pakai QDA)
Perbedaan vektor rata-rata (MANOVA)	p < 0,0001	Kelompok berbeda nyata

Pembangunan Model LDA dan QDA

Model LDA

lda(Species ~ ., data = latih) membangun fungsi diskriminan linear dengan Species sebagai kelompok dan semua peubah lain (tanda .) sebagai prediktor. Model dibangun hanya dari data latih. Setelah itu sd_pooled() dipakai untuk mendapatkan koefisien terstandardisasi (koefisien asli dikali simpangan baku gabungan), yang memungkinkan membandingkan besar kontribusi tiap peubah secara adil walau satuannya berbeda.

lda_fit <- lda(Species ~ ., data = latih)   
lda_fit

## Call:
## lda(Species ~ ., data = latih)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa         5.020000    3.397143     1.477143   0.2457143
## versicolor     5.988571    2.768571     4.271429   1.3171429
## virginica      6.651429    2.994286     5.608571   2.0600000
## 
## Coefficients of linear discriminants:
##                     LD1        LD2
## Sepal.Length  0.7368248  0.1726288
## Sepal.Width   1.2334070 -2.2628483
## Petal.Length -2.1311089  0.7858110
## Petal.Width  -2.5226526 -2.6041711
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9894 0.0106

# Koefisien terstandardisasi: koefisien asli x simpangan baku gabungan
koef_std <- lda_fit$scaling * sd_pooled(X_latih, g_latih)
kable(round(koef_std, 3), caption = "Koefisien fungsi diskriminan terstandardisasi (LDA)")

Koefisien fungsi diskriminan terstandardisasi (LDA)
	LD1	LD2
Sepal.Length	0.358	0.084
Sepal.Width	0.423	-0.776
Petal.Length	-0.903	0.333
Petal.Width	-0.540	-0.558

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

Prior probability ketiga spesies sama (0,333), karena data latih seimbang.
Group means menunjukkan pola yang sama dengan eksplorasi: nilai petal naik dari setosa ke versicolor ke virginica.
Karena ada 3 kelompok, terbentuk 2 fungsi diskriminan (jumlah kelompok − 1). LD1 menjelaskan 98,94% pemisahan antarkelompok, sedangkan LD2 hanya 1,06%. Jadi hampir seluruh pembeda spesies tercakup dalam LD1.
Pada koefisien terstandardisasi LD1, peubah dengan nilai mutlak terbesar adalah Petal.Length (-0,903). Peubah petal mendominasi, sejalan dengan hasil eksplorasi bahwa petal adalah pembeda terkuat.

skor_latih <- data.frame(predict(lda_fit)$x, Species = latih$Species)  

ggplot(skor_latih, aes(x = LD1, y = LD2, colour = Species)) +
  geom_point(size = 2.2, alpha = 0.85) +
  stat_ellipse(level = 0.95, linewidth = 0.7) +
  scale_colour_manual(values = warna) +
  labs(title = "Data latih pada ruang diskriminan (LD1 vs LD2)") +
  theme_minimal(base_size = 11)

Sumbu LD1 (horizontal) memisahkan setosa jauh di salah satu ujung, sedangkan versicolor dan virginica terpisah lebih dekat tetapi masih terbedakan dengan baik. Gerombol yang hampir tidak bersinggungan menandakan fungsi diskriminan berhasil memisahkan spesies; titik-titik pada daerah singgung versicolor-virginica adalah kandidat salah klasifikasi.

Model QDA

qda() dipanggil dengan cara yang sama seperti lda(). Bedanya, QDA menaksir matriks kovarians terpisah untuk tiap spesies, sehingga batas keputusannya berupa kurva dan tidak memiliki koefisien tunggal seperti LDA.

qda_fit <- qda(Species ~ ., data = latih)   
qda_fit

## Call:
## qda(Species ~ ., data = latih)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa         5.020000    3.397143     1.477143   0.2457143
## versicolor     5.988571    2.768571     4.271429   1.3171429
## virginica      6.651429    2.994286     5.608571   2.0600000

Prior dan rata-rata kelompok sama dengan LDA. Perbedaannya tidak tampak pada output ini, melainkan pada cara model memperlakukan sebaran tiap spesies: QDA mengizinkan setiap spesies memiliki bentuk dan ukuran sebaran sendiri. Karena itu kualitas QDA dinilai lewat hasil klasifikasi, bukan lewat koefisien.

Evaluasi Model pada Data Latih

predict(model) tanpa newdata mengklasifikasikan ulang data latih. Elemen $class berisi spesies prediksi. Prediksi dibandingkan dengan label asli: mean(prediksi != aktual) adalah tingkat kesalahan (APER), mean(prediksi == aktual) adalah akurasi, dan diag(m) / rowSums(m) pada matriks klasifikasi memberi sensitivitas tiap spesies.

pred_lda_latih <- predict(lda_fit)$class     
pred_qda_latih <- predict(qda_fit)$class     

err_latih <- data.frame(
  Model             = c("LDA", "QDA"),
  `Jumlah salah`    = c(sum(pred_lda_latih != g_latih), sum(pred_qda_latih != g_latih)),
  `Tingkat kesalahan (APER)` = c(mean(pred_lda_latih != g_latih), mean(pred_qda_latih != g_latih)),
  Akurasi           = c(mean(pred_lda_latih == g_latih), mean(pred_qda_latih == g_latih)),
  check.names = FALSE
)
kable_d(err_latih, 3:4, 4, caption = "Kemampuan klasifikasi pada data latih")

Kemampuan klasifikasi pada data latih
Model	Jumlah salah	Tingkat kesalahan (APER)	Akurasi
LDA	2	0,0190	0,9810
QDA	2	0,0190	0,9810

kable(tabel_konfusi(g_latih, pred_lda_latih), caption = "Matriks klasifikasi LDA (data latih)")

Matriks klasifikasi LDA (data latih)
Aktual ↓ / Prediksi →	setosa	versicolor	virginica
setosa	35	0	0
versicolor	0	34	1
virginica	0	1	34

kable(tabel_konfusi(g_latih, pred_qda_latih), caption = "Matriks klasifikasi QDA (data latih)")

Matriks klasifikasi QDA (data latih)
Aktual ↓ / Prediksi →	setosa	versicolor	virginica
setosa	35	0	0
versicolor	0	34	1
virginica	0	1	34

# Sensitivitas: persentase tiap spesies yang dikenali dengan benar
sens_latih <- data.frame(
  Spesies = levels(g_latih),
  LDA = diag(table(g_latih, pred_lda_latih)) / table(g_latih),
  QDA = diag(table(g_latih, pred_qda_latih)) / table(g_latih)
)[, c(1, 3, 5)]
names(sens_latih) <- c("Spesies", "Sensitivitas LDA", "Sensitivitas QDA")
kable_d(sens_latih, 2:3, 3, row.names = FALSE, caption = "Sensitivitas per spesies (data latih)")

Sensitivitas per spesies (data latih)
Spesies	Sensitivitas LDA	Sensitivitas QDA
setosa	1,000	1,000
versicolor	0,971	0,971
virginica	0,971	0,971

Pada data latih, LDA salah mengklasifikasikan 2 dari 105 amatan (APER 0,0190) dan QDA 2 dari 105 (APER 0,0190). Akurasi keduanya 98,1% dan 98,1%. Pada matriks klasifikasi, setosa dikenali sempurna, sedangkan kesalahan yang ada hanya terjadi antara versicolor dan virginica, dua spesies yang ukurannya paling berdekatan. Ini sejalan dengan grafik sebelumnya.

Perlu diingat bahwa angka pada data latih bersifat optimistis, karena data yang sama dipakai untuk membangun sekaligus menguji model. Penilaian yang lebih jujur ada pada data uji.

Evaluasi Model pada Data Uji

Sekarang model yang sudah dibangun dipakai untuk memprediksi data uji lewat argumen newdata = uji. Selain $class, predict() menyediakan $posterior, yaitu peluang posterior tiap spesies untuk setiap amatan. Peluang ini dipakai untuk melihat seberapa yakin model pada amatan yang salah klasifikasi.

g_uji <- uji$Species                                   

pred_lda_uji <- predict(lda_fit, newdata = uji)        
pred_qda_uji <- predict(qda_fit, newdata = uji)        
kelas_lda_uji <- pred_lda_uji$class
kelas_qda_uji <- pred_qda_uji$class

err_uji <- data.frame(
  Model             = c("LDA", "QDA"),
  `Jumlah salah`    = c(sum(kelas_lda_uji != g_uji), sum(kelas_qda_uji != g_uji)),
  `Tingkat kesalahan` = c(mean(kelas_lda_uji != g_uji), mean(kelas_qda_uji != g_uji)),
  Akurasi           = c(mean(kelas_lda_uji == g_uji), mean(kelas_qda_uji == g_uji)),
  check.names = FALSE
)
kable_d(err_uji, 3:4, 4, caption = "Kemampuan klasifikasi pada data uji")

Kemampuan klasifikasi pada data uji
Model	Jumlah salah	Tingkat kesalahan	Akurasi
LDA	1	0,0222	0,9778
QDA	1	0,0222	0,9778

kable(tabel_konfusi(g_uji, kelas_lda_uji), caption = "Matriks klasifikasi LDA (data uji)")

Matriks klasifikasi LDA (data uji)
Aktual ↓ / Prediksi →	setosa	versicolor	virginica
setosa	15	0	0
versicolor	0	14	1
virginica	0	0	15

kable(tabel_konfusi(g_uji, kelas_qda_uji), caption = "Matriks klasifikasi QDA (data uji)")

Matriks klasifikasi QDA (data uji)
Aktual ↓ / Prediksi →	setosa	versicolor	virginica
setosa	15	0	0
versicolor	0	14	1
virginica	0	0	15

# Amatan uji yang salah diklasifikasikan oleh LDA dan/atau QDA
idx_salah <- which(kelas_lda_uji != g_uji | kelas_qda_uji != g_uji)

if (length(idx_salah) > 0) {
  p_lda <- apply(pred_lda_uji$posterior[idx_salah, , drop = FALSE], 1, max)  
  p_qda <- apply(pred_qda_uji$posterior[idx_salah, , drop = FALSE], 1, max)  
  kable(data.frame(
    `Baris asli` = rownames(uji)[idx_salah],
    Aktual       = g_uji[idx_salah],
    `Prediksi LDA (peluang)` = paste0(kelas_lda_uji[idx_salah], " (", id(p_lda, 2), ")"),
    `Prediksi QDA (peluang)` = paste0(kelas_qda_uji[idx_salah], " (", id(p_qda, 2), ")"),
    check.names = FALSE, row.names = NULL),
    caption = "Amatan uji yang salah diklasifikasikan")
} else {
  cat("Tidak ada amatan uji yang salah diklasifikasikan.\n")
}

Amatan uji yang salah diklasifikasikan
Baris asli	Aktual	Prediksi LDA (peluang)	Prediksi QDA (peluang)
71	versicolor	virginica (0,76)	virginica (0,87)

skor_uji <- data.frame(predict(lda_fit, newdata = uji)$x, Species = g_uji)  
skor_uji$Salah <- kelas_lda_uji != g_uji | kelas_qda_uji != g_uji           

ggplot(skor_uji, aes(x = LD1, y = LD2, colour = Species)) +
  geom_point(size = 2.6) +
  geom_point(data = subset(skor_uji, Salah), shape = 21, size = 6, stroke = 1.1,
             colour = "black", fill = NA) +                         
  scale_colour_manual(values = warna) +
  labs(title = "Data uji pada ruang diskriminan (LD1 vs LD2)") +
  theme_minimal(base_size = 11)

Pada data uji (45 amatan), LDA salah pada 1 amatan (tingkat kesalahan 0,0222, akurasi 97,78%) dan QDA salah pada 1 amatan (tingkat kesalahan 0,0222, akurasi 97,78%). Seluruh setosa terklasifikasi benar. Amatan yang salah terletak pada daerah singgung versicolor dan virginica pada grafik, dan peluang posterior pada tabel kesalahan jauh dari 1, yang menandakan model memang ragu pada amatan batas tersebut. Karena jumlah kesalahan LDA dan QDA sama, perbedaan kinerja kedua model tidak dapat dibedakan pada data uji ini.

Kestabilan Hasil

Hasil di atas berasal dari satu kali pembagian data (seed 2026), sehingga bisa saja kebetulan. Dua pemeriksaan tambahan dilakukan:

Pembagian berulang. Fungsi satu_pembagian() mengulang seluruh proses (bagi 70:30 berstrata, bangun LDA dan QDA, hitung kesalahan uji). replicate(200, ...) menjalankannya 200 kali.
Validasi silang LOOCV (leave-one-out). Dengan CV = TRUE, setiap bunga diklasifikasikan oleh model yang dibangun dari 149 bunga lainnya, sehingga seluruh 150 data dipakai untuk menguji tanpa bocoran.

satu_pembagian <- function() {
  idx <- sort(unlist(lapply(split(seq_len(nrow(iris_data)), iris_data$Species),
                            function(ii) sample(ii, size = round(0.7 * length(ii))))))
  tr <- iris_data[idx, ];  te <- iris_data[-idx, ]
  c(LDA = mean(predict(lda(Species ~ ., tr), te)$class != te$Species),   
    QDA = mean(predict(qda(Species ~ ., tr), te)$class != te$Species))   
}

set.seed(1)                                            
ulang <- as.data.frame(t(replicate(200, satu_pembagian())))

ringkas_ulang <- data.frame(
  Model = c("LDA", "QDA"),
  `Rata-rata` = c(mean(ulang$LDA), mean(ulang$QDA)),
  `Simpangan baku` = c(sd(ulang$LDA), sd(ulang$QDA)),
  Minimum = c(min(ulang$LDA), min(ulang$QDA)),
  Maksimum = c(max(ulang$LDA), max(ulang$QDA)),
  check.names = FALSE
)
kable_d(ringkas_ulang, 2:5, 4, caption = "Tingkat kesalahan data uji pada 200 pembagian acak")

Tingkat kesalahan data uji pada 200 pembagian acak
Model	Rata-rata	Simpangan baku	Minimum	Maksimum
LDA	0,0211	0,0166	0,0000	0,0667
QDA	0,0248	0,0193	0,0000	0,0667

# Berapa kali QDA lebih baik, sama, atau lebih buruk dari LDA?
kable(data.frame(
  Hasil = c("QDA lebih baik", "Sama", "LDA lebih baik"),
  `Jumlah dari 200 pembagian` = c(sum(ulang$QDA < ulang$LDA), sum(ulang$QDA == ulang$LDA),
                                  sum(ulang$QDA > ulang$LDA)),
  check.names = FALSE), caption = "Perbandingan LDA dan QDA pada tiap pembagian")

Perbandingan LDA dan QDA pada tiap pembagian
Hasil	Jumlah dari 200 pembagian
QDA lebih baik	18
Sama	133
LDA lebih baik	49

# LOOCV pada seluruh 150 data
cv_lda <- lda(Species ~ ., data = iris_data, CV = TRUE)$class
cv_qda <- qda(Species ~ ., data = iris_data, CV = TRUE)$class
kable_d(data.frame(Model = c("LDA", "QDA"),
                   `Tingkat kesalahan LOOCV` = c(mean(cv_lda != iris_data$Species),
                                                 mean(cv_qda != iris_data$Species)),
                   check.names = FALSE), 2, 4, caption = "Tingkat kesalahan LOOCV (150 data)")

Tingkat kesalahan LOOCV (150 data)
Model	Tingkat kesalahan LOOCV
LDA	0,0200
QDA	0,0267

ulang %>%
  pivot_longer(everything(), names_to = "Model", values_to = "Kesalahan") %>%
  ggplot(aes(x = Model, y = Kesalahan, fill = Model)) +
  geom_boxplot(alpha = 0.6, show.legend = FALSE, outlier.shape = NA) +
  geom_jitter(width = 0.12, height = 0, alpha = 0.25, size = 1, show.legend = FALSE) +
  scale_fill_manual(values = c(LDA = "#3182ce", QDA = "#dd6b20")) +
  labs(title = "Kestabilan tingkat kesalahan data uji", x = NULL, y = "Tingkat kesalahan") +
  theme_minimal(base_size = 11)

Dari 200 pembagian acak, rata-rata tingkat kesalahan data uji adalah 0,0211 untuk LDA dan 0,0248 untuk QDA. QDA lebih baik pada 18 pembagian, sama pada 133, dan LDA lebih baik pada 49. Pada LOOCV, kesalahan LDA 0,0200 dan QDA 0,0267.

Kedua pemeriksaan sama-sama menunjukkan LDA sedikit lebih baik daripada QDA, walaupun uji Box’s M menolak kesamaan ragam. Kemungkinan penyebabnya: QDA harus menaksir matriks kovarians terpisah untuk tiap spesies (lebih banyak parameter) dari hanya 35 amatan per spesies, sehingga taksirannya lebih tidak stabil. Pada data kecil, model yang lebih sederhana sering menang walaupun asumsinya tidak sempurna. Yang terpenting, kedua model sama-sama sangat akurat (kesalahan hanya beberapa persen) dan selisihnya kecil, sehingga angka dari satu kali pembagian layak dipercaya sebagai gambaran umum.

Kesimpulan

tabel_ringkas_akhir <- data.frame(
  Model = c("LDA", "QDA"),
  `Kesalahan latih` = err_latih[, 3],
  `Kesalahan uji`   = err_uji[, 3],
  `Akurasi uji`     = err_uji[, 4],
  `Kesalahan LOOCV` = c(mean(cv_lda != iris_data$Species), mean(cv_qda != iris_data$Species)),
  `Rata-rata 200 pembagian` = c(mean(ulang$LDA), mean(ulang$QDA)),
  check.names = FALSE)
kable_d(tabel_ringkas_akhir, 2:6, 4, caption = "Ringkasan kinerja LDA dan QDA")

Ringkasan kinerja LDA dan QDA
Model	Kesalahan latih	Kesalahan uji	Akurasi uji	Kesalahan LOOCV	Rata-rata 200 pembagian
LDA	0,0190	0,0222	0,9778	0,0200	0,0211
QDA	0,0190	0,0222	0,9778	0,0267	0,0248

Asumsi. Normalitas multivariat tidak ditolak pada ketiga spesies (Mardia), sedangkan matriks ragam-peragam antarspesies tidak sama (Box’s M, p < 0,0001). Dengan demikian, secara teori QDA lebih sesuai. MANOVA menunjukkan rata-rata ketiga spesies berbeda nyata, sehingga analisis diskriminan layak dilakukan.
Pembeda utama. Panjang dan lebar petal adalah peubah yang paling membedakan spesies; LD1 menjelaskan 98,9% pemisahan antarkelompok.
Kinerja. Pada data latih, LDA dan QDA salah 2 dan 2 dari 105 amatan. Pada data uji, keduanya salah 1 dan 1 dari 45 amatan (akurasi 97,78% dan 97,78%). Kesalahan hanya terjadi antara versicolor dan virginica; setosa selalu terklasifikasi benar.
Pemilihan model. Kedua model memberi kinerja yang setara pada data uji, dan selisihnya tetap kecil pada pengulangan pembagian maupun LOOCV. Secara teori QDA lebih sesuai (kesamaan ragam ditolak), tetapi secara empiris LDA sedikit lebih stabil pada data sekecil ini. Maka LDA dapat dipakai sebagai model praktis, dengan QDA sebagai pembanding; keduanya memberi hasil klasifikasi yang hampir sama.

Keterbatasan

Data uji hanya 45 amatan, sehingga satu amatan salah saja mengubah tingkat kesalahan sekitar 2,2 poin persen; estimasi kesalahan uji memiliki ketidakpastian yang cukup besar.
Pembagian 70:30 pada analisis utama hanya dilakukan sekali sehingga bergantung pada seed. Hal ini sudah diperiksa pada bagian Kestabilan Hasil.
Evaluasi pada data latih bersifat optimistis dan tidak dijadikan dasar utama penilaian.
Dengan 35 amatan per spesies pada data latih, daya uji Mardia dan Box’s M terbatas, dan Box’s M peka terhadap pelanggaran normalitas.
Dataset iris adalah data yang relatif “mudah” (kelompok terpisah baik), sehingga performa tinggi di sini belum tentu berlaku pada data lain.

Analisis Diskriminan pada Dataset Iris

Tugas Praktikum Teknik Peubah Ganda

Fadly Syabani

09-11-2026