The EM algorithm by example of k-means clustering

Mungkin algoritma yang paling terkenal untuk mengelompokkan pengamatan ke dalam kelompok adalah algoritma k-means. Kita akan melihat bahwa algoritma ini hanyalah sebuah varian dari algoritma EM. Diberikan n objek, yang dicirikan oleh p variabel, kita ingin mempartisi mereka ke dalam nC cluster {C1,C2,…,CnC} sedemikian rupa sehingga cluster Ck memiliki n(k) anggota dan setiap observasi berada dalam satu cluster. Vektor rata-rata (pusat, prototipe), Vk , dari sebuah klaster Ck didefinisikan sebagai pusat (centroid) dari klaster tersebut dan komponen-komponen dari vektor rata-rata tersebut dapat dihitung dengan:

vk(∈Rp)=(1n(k)∑i=1n(k)x(k)i1,…,1n(k)∑i=1n(k)x(k)ip)

Dimana n(k) adalah jumlah pengamatan dalam cluster Ck dan x(k)i adalah pengamatan ke-i yang termasuk dalam cluster Ck . Untuk setiap cluster C1,…,CnC , cluster yang sesuai berarti V={v1,…,vnC} dihitung.

Kita juga perlu menentukan jumlah cluster dalam partisi keluaran. Dimulai dari lokasi awal yang diberikan dari centroid cluster nC , algoritme menggunakan titik-titik data untuk secara iteratif merelokasi centroid dan mengalokasikan kembali titik-titik ke centroid terdekat. Proses ini terdiri dari langkah-langkah berikut:

Pilih partisi awal dengan cluster nC . Langkah 3. Langkah E: (Hitung ulang pusat-pusat klaster menggunakan keanggotaan klaster saat ini. Langkah 3. Langkah M: Tetapkan setiap objek ke pusat klaster terdekat, dengan keanggotaan baru. Lanjutkan ke langkah 2 hingga keanggotaan klaster dan pusat klaster tidak berubah di luar batas yang ditentukan. K-means clustering mengoptimalkan fungsi objektif: Translated with DeepL.com (free version)

J(X,V,U)=∑k=1nc∑i=1nuikd2(xi,vk) Dimana X={x1,…,xn} adalah kumpulan data dengan observasi dan variabel, v1,…,vnc adalah matriks pusat klaster (prototipe) berdimensi nc . U=[(uik)] adalah matriks dengan koefisien keanggotaan uik untuk pengamatan xi pada sebuah klaster nc . Oleh karena itu, U berdimensi n . d adalah jarak Euclidean antara pengamatan dan pusat cluster. nc menentukan jumlah cluster.

Algoritma k-means dapat diimplementasikan sebagai berikut. Tetapkan nc , 2≤nc<n , dan pilih toleransi penghentian δ>0 , sebagai contoh, 0.001 . Inisialisasi U(0) (misalnya, secara acak).

REPEAT untuk $ r = 1, 2, titik-titik$

LANGKAH E: Hitung pusat-pusat gugus: vk=1∑ni=1uik(∑i=1nuikxi1,…,∑i=1nuikxip),for k=1,…,nc 2. M-langkah: Perbarui U(r) : Alokasikan kembali keanggotaan klaster.

u(r)ij={10if d(xi,v(r)j)=min1≤l≤ncd(xi,v(r)l)lainnya

UNTIL $ | ^{(r)} - ^{(r-1)} | < $

Kami mendefinisikan seluruh algoritma k-means sebagai berikut. Ini hanya untuk memahami kode dan algoritmanya. Untuk implementasi profesional dari k-means, lihat misalnya, (Leisch 2006).

Untuk algoritma cluster kita membutuhkan fungsi jarak. Kami menggunakan Jarak Manhattan:

distMan <- function(x, centers){
if(class(x) == "data.frame") x <- as.matrix(x)
d <- matrix(0, nrow=nrow(x), ncol=nrow(centers))
## dist center to observations for each cluster
for(k in 1:nrow(centers)){
d[,k] <- abs( colSums((t(x) - centers[k,])) )
}
return(d)
}

Dan kita membutuhkan fungsi yang menghitung rata-rata, misalnya, kita dapat menggunakan aritmatika rata-rata, tetapi kita juga dapat menggunakan median seperti yang berikut ini:

means <- function(x, cluster){
cen <- matrix(NA, nrow=max(cluster), ncol <- ncol(x))
## cluster means for each cluster
for(n in 1:max(cluster)){
cen[n,] <- apply(x[cluster==n,], 2, median)
}
return(cen)
}

Kita menulis sebuah fungsi untuk algoritma k-means, yang mengimplementasikan rumus-rumus sebelumnya. Untuk membuat beberapa plot yang menunjukkan pendekatan EM secara detail, kita melakukannya dalam sebuah perulangan for :

my_kmeans <- function(x, k, clmeans = means, distance = distMan, iter =
99999, seed = NULL){
if(!is.null(seed)) set.seed(seed)
cent <- newcent <- x[sample(1:nrow(x), size=k), ]
oldclust <- 0
j <- 0
for(i in 1:iter){ # better: while()
j <- j + 1
cent <- newcent
## M-step
dist <- distance(x, cent)
clust <- max.col(-dist)
## E-step
newcent <- clmeans(x, clust)
if(all(clust == oldclust)) break()
oldclust <- clust
}
res <- list(centers = cent,
cluster = factor(clust),
iterations = j)
return(res)
}

Seperti yang bisa kita lihat, dalam pengelompokan k-means, langkah E adalah langkah penyesuaian dan langkah M adalah langkah penugasan. Mengulangi langkah E dan M secara berulang akan meningkatkan solusi. Ini berarti bahwa J(X,V,U) menjadi lebih kecil di setiap iterasi. Kita hentikan algoritmanya jika penugasan klaster tidak berubah lagi. Mari kita ambil beberapa data. Kita ingin membuatnya tetap sederhana dan kita ingin menampilkan klaster-klaster secara visual. Jadi kita telah mengambil satu set data dua dimensi dan menunjukkannya pada Gambar 9.1.

data(Nclus, package = "flexclust")
x <- data.frame(Nclus)
library("ggplot2")
qplot(X1, X2, data=data.frame(Nclus))
## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Berikut ini kami memplotkan hasil setelah iterasi 1, 2, dan setelah konvergensi. Alih-alih implementasi algoritma k-means yang disederhanakan, kami menggunakan implementasi k-means default dari R. Beberapa varian dari k-means ada, di mana kami memilih algoritma “MacQueen”, tetapi hanya untuk mengeksplorasi algoritma tersebut (metode default, “Hartigan-Wong” konvergen terlalu cepat untuk menunjukkan langkah-langkah algoritma). Perhatikan bahwa algoritma k-means dimulai dengan pusat-pusat klaster yang dipilih secara acak. Oleh karena itu, kita harus mengatur seed untuk memastikan permulaan yang sama di setiap pemanggilan k-means:

set.seed(123456)
cl1 <- kmeans(x, centers = 4, iter.max = 1, algorithm = "MacQueen")
## Warning: did not converge in 1 iteration
set.seed(123456)
cl2 <- kmeans(x, centers = 4, iter.max = 2, algorithm = "MacQueen")
## Warning: did not converge in 2 iterations
set.seed(123456)
cl3 <- kmeans(x, centers = 4, iter.max = 3, algorithm = "MacQueen")
## Warning: did not converge in 3 iterations
set.seed(123456)
cl4 <- kmeans(x, centers = 4, algorithm = "MacQueen")


cl1$centers
##           X1         X2
## 1  1.4080304  1.9778530
## 2 -1.9756156  6.0385730
## 3  0.1149815 -0.7912543
## 4  5.7943888  2.6065412

Perhatikan bahwa k-means hanya memperhitungkan pusat-pusatnya dan bekerja dengan fungsi jarak untuk menghitung jarak dari pengamatan ke pusat-pusat cluster. Pendekatan lain (yang lebih baik) adalah dengan memasukkannya ke dalam bentuk cluster. Hal ini diimplementasikan dalam kerangka kerja pengelompokan berbasis model (Fraley dan Raftery 2002). Prosedur berbasis model sebagian besar memberikan hasil pengelompokan yang lebih baik (Templ, Filzmoser, dan Reimann 2008), namun secara komputasi lebih rumit karena pada setiap langkah E, kovariansi dari setiap klaster harus dihitung.

The EM algorithm for the imputation of missing values

Algoritma EM digunakan secara luas untuk imputasi nilai yang hilang. Implementasinya meliputi (van Buuren dan Groothuis-Oudshoorn 2011), (Schafer 1997), (Templ, Alfons, dan Filzmoser 2011), (Raghunathan dkk. 2001), dan (Gelman dan Hill 2011). Berikut ini kami ingin menunjukkan bagaimana algoritma EM bekerja secara umum untuk masalah-masalah seperti ini. Pertama, kami mengambil satu set data untuk diimputasi. Kami memilih lagi data tidur:

library("MASS")
library("robustbase")
library("VIM")
## Loading required package: colorspace
## Loading required package: grid
## VIM is ready to use.
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
## 
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
## 
##     sleep
data("sleep")
str(sleep)
## 'data.frame':    62 obs. of  10 variables:
##  $ BodyWgt : num  6654 1 3.38 0.92 2547 ...
##  $ BrainWgt: num  5712 6.6 44.5 5.7 4603 ...
##  $ NonD    : num  NA 6.3 NA NA 2.1 9.1 15.8 5.2 10.9 8.3 ...
##  $ Dream   : num  NA 2 NA NA 1.8 0.7 3.9 1 3.6 1.4 ...
##  $ Sleep   : num  3.3 8.3 12.5 16.5 3.9 9.8 19.7 6.2 14.5 9.7 ...
##  $ Span    : num  38.6 4.5 14 NA 69 27 19 30.4 28 50 ...
##  $ Gest    : num  645 42 60 25 624 180 35 392 63 230 ...
##  $ Pred    : int  3 3 1 5 3 4 1 4 1 1 ...
##  $ Exp     : int  5 1 1 2 5 4 1 5 2 1 ...
##  $ Danger  : int  3 3 1 3 4 4 1 4 1 1 ...

Missing values are included in some of the all variables, like in variable Sleep:

apply(sleep, 2, function(x) any(is.na(x)))
##  BodyWgt BrainWgt     NonD    Dream    Sleep     Span     Gest     Pred 
##    FALSE    FALSE     TRUE     TRUE     TRUE     TRUE     TRUE    FALSE 
##      Exp   Danger 
##    FALSE    FALSE

Bagaimana kita dapat mengimputasi nilai yang hilang dalam variabel Tidur? Hal ini dapat dilakukan dengan melakukan pencocokan regresi. Idealnya, kita sudah menginisialisasi nilai yang hilang dari algoritme yang berkinerja baik untuk imputasi, misalnya, pendekatan imputasi tetangga terdekat (Templ, Alfons, dan Filzmoser 2011). Namun, untuk melihat kemajuan dari EM, kita mulai dengan inisialisasi nilai yang hilang di variabel Sleep dengan sangat buruk, menggunakan inisialisasi terburuk dengan jumlah yang besar. Perhatikan bahwa kita juga membutuhkan indeks dari nilai-nilai yang hilang:

## index of missing values
ind <- data.frame(is.na(sleep))
## initialization
sleep <- kNN(sleep)
## Time difference of 0.04911399 secs
## overwrite missing initialization with bad choice
sleep$Sleep[ind$Sleep] <- 2240 # bad initialization
## initialized missing values in variable sleep
sleep$Sleep[ind$Sleep]
## [1] 2240 2240 2240 2240

We then fit a model for the first variable. The model results (regression coefficients) are then used to predict the missing values

## E-step (1)
lm1 <- lm(Sleep ~ log(BodyWgt) + log(BrainWgt) + NonD + Danger, data =
sleep)
## M-step (1)
sleep$Sleep[ind$Sleep] <- predict(lm1)[ind$Sleep]
## print of updated missing values
sleep$Sleep[ind$Sleep]
## [1] 469.5127 559.9771 408.6845 229.0985

Kita dapat melanjutkan dengan iterasi kedua.

## E-step (2)
lm1 <- lm(Sleep ~ log(BodyWgt) + log(BrainWgt) + NonD + Danger, data =
sleep)
## M-step (2)
sleep$Sleep[ind$Sleep] <- predict(lm1)[ind$Sleep]
## print of updated missing values
sleep$Sleep[ind$Sleep]
## [1] 101.9265 121.6146  90.1618  48.7181

Kita melihat bahwa nilainya masih banyak berubah. Mari kita lakukan ini secara iteratif sampai nilainya tidak akan berubah lebih dari ambang batas yang sangat kecil. Kita bisa menulis sebuah fungsi kecil untuk imputasi variabel Sleep:

EMregSleep <- function(method = lm, eps = 0.001, init = "bad"){
## index of missing values
ind <- is.na(sleep)
colnames(ind) <- colnames(sleep)
indsleep <- ind[, "Sleep"]
## initialization
if(init == "bad"){
sleep <- kNN(sleep, imp_var = FALSE)
sleep$Sleep[indsleep] <- 2240 # bad initialization
}
if(init == "worst"){
sleep[ind] <- 2240 # worst initialization
}
iteration <- 0
criteria <- 99999
while(criteria > eps){
iteration <- iteration + 1
prev_sol <- sleep$Sleep[indsleep]
## E-step
lm1 <- method(Sleep ~ log(BodyWgt) + log(BrainWgt) + NonD + Danger,
data = sleep)
## M-step
sleep$Sleep[indsleep] <- predict(lm1)[indsleep]
criteria <- sqrt(sum((prev_sol - sleep$Sleep[indsleep])^2))
}
res <- list("imputed" = sleep,
"iteration" = iteration,
lastmodel = lm1)
return(res)
}

Sekali lagi kita memuat set data tidur, memperhitungkannya dan melihat nilai yang diperhitungkan dalam variabel Tidur:

data("sleep")
sleepImp <- EMregSleep()
## Time difference of 0.179677 secs
missVals <- sleepImp$imputed$Sleep[ind$Sleep]
missVals
## [1]  3.845778 13.122764  3.658173 16.975766
sleepImp$iteration
## [1] 11

Namun, kami memperhitungkan dengan nilai yang diharapkan, yang akan menyebabkan penurunan varians karena kami tidak memperhitungkan ketidakpastian dan distribusi nilai yang hilang (bandingkan estimasi varians dengan nilai yang hilang di Bab 8, Aplikasi Metode Resampling dan Uji Monte Carlo). Untuk mempertimbangkan varians, kami mengambil sampel residual (bandingkan pendekatan regresi residual bootstrapping di Bab 8, Aplikasi Metode Resampling dan Uji Monte Carlo):

missVals + sample(residuals(sleepImp$lastmodel), length(missVals))
##         7        45        54        36 
##  4.122294 13.199528  3.844590 17.827567

Sebelumnya kita telah melihat bahwa kita membutuhkan 11 iterasi. Selain itu, model regresi OLS mungkin juga dipengaruhi oleh outlier, sehingga lebih baik menggantinya dengan metode arobust. Kita sudah melihat hasil yang baik setelah iterasi pertama (tidak ditampilkan di sini). Kami mendapatkan hasil yang sedikit berbeda yang biasanya lebih dapat dipercaya daripada menggunakan metode non-robust (Templ, Kowarik, dan Filzmoser 2011). Hasil OLS dapat menjadi rusak, terutama dengan adanya pencilan dalam prediktor. Namun, kami melihat bahwa bahkan dengan inisialisasi terburuk (juga pencilan yang sangat besar pada prediktor) hasilnya tetap terlihat baik (meskipun kami lebih memilih metode robust):

data("sleep")
## OLS regression
lm_ols <- EMregSleep(method = lm, init = "worst")
## M-estimation
lm_rlm <- EMregSleep(method = rlm, init = "worst", eps= 0.01)
lm_ols$imputed$Sleep[ind[, "Sleep"]]
## [1]  4.239191  8.169014  4.368256 13.775087
## [1] 4.239191 8.169014 4.368256 13.775087
lm_rlm$imputed$Sleep[ind[, "Sleep"]]
## [1]  3.766792  7.788943  3.925772 13.700029

Dari angka-angka ini kita dapat melihat bahwa hasil OLS sangat dipengaruhi oleh outlier. Dibandingkan dengan estimasi sebelumnya (menggunakan inisialisasi yang buruk, bukan yang terburuk), nilai yang diimputasi terlalu tinggi. Hal ini tidak seekstrim ketika menggunakan estimator-M, tetapi dibandingkan dengan implementasi dalam fungsi irmi (lihat contoh berikut), kami meremehkan nilai kedua dan keempat.

Kita telah membahas bagaimana cara mengimputasi satu variabel, tetapi secara umum kita ingin mengimputasi semua variabel dalam kumpulan data. Kumpulan data juga dapat terdiri dari tidak hanya variabel kontinu tetapi juga campuran variabel kontinu, semi-kontinu, kategorik, biner, dan/atau cacahan. Imputasi berbasis EM yang kuat memperhitungkan hal ini (dan lebih banyak lagi, seperti menentukan model untuk setiap variabel) dan diimplementasikan dalam fungsi irmi (Templ, Kowarik, dan Filzmoser 2011) di dalam paket R VIM (Templ, Alfons, dan Filzmoser 2011):

data("sleep")
sleepImp <- irmi(sleep)
## Time difference of 0.03798294 secs
sleepImp[ind[, "Sleep"], "Sleep"]
## [1]  3.748899 10.089591  3.156300 17.085060

Kami melihat ini sangat dekat dengan solusi awal di mana kami mengambil inisialisasi yang lebih baik dari nilai yang hilang. Ini merupakan indikasi keberhasilan irmi. Kita dapat menggunakan metode lain, seperti mice (karena irmi biasanya digunakan untuk imputasi berganda):

library("mice")
## 
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
## 
##     filter
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
## Loading required package: Rcpp
## mice 2.25 2015-11-09
data("sleep")
em_mice <- mice(sleep, m = 1)
## 
##  iter imp variable
##   1   1  NonD  Dream  Sleep  Span  Gest
##   2   1  NonD  Dream  Sleep  Span  Gest
##   3   1  NonD  Dream  Sleep  Span  Gest
##   4   1  NonD  Dream  Sleep  Span  Gest
##   5   1  NonD  Dream  Sleep  Span  Gest
## Warning: Number of logged events: 2
em_mice$imp$Sleep
##       1
## 21 12.8
## 31  9.8
## 41  2.6
## 62 11.2
## now with bad intitialisation in predictors
sleep[is.na(sleep)] <- 2240
sleep$Sleep[ind[, "Sleep"]] <- NA
em_mice <- mice(sleep, m = 1)
## 
##  iter imp variable
##   1   1  Sleep
##   2   1  Sleep
##   3   1  Sleep
##   4   1  Sleep
##   5   1  Sleep
## Warning: Number of logged events: 5
em_mice$imp$Sleep
##       1
## 21  3.9
## 31  2.6
## 41  3.1
## 62 13.8

Kami melihat bahwa kami mendapatkan hasil yang sama sekali berbeda, segera setelah pencilan hadir dalam kumpulan data. Hal ini juga berlaku untuk implementasi algoritma EM lainnya, yang hanya cocok jika set data mendekati normal multivariat. Segera setelah hal ini dilanggar (seperti yang biasanya terjadi dalam praktiknya), irmi mungkin merupakan pilihan yang baik. Perhatikan bahwa dalam irmi banyak parameter lain yang dapat ditentukan (tidak dibahas di sini):

args(irmi)
## function (x, eps = 5, maxit = 100, mixed = NULL, mixed.constant = NULL, 
##     count = NULL, step = FALSE, robust = FALSE, takeAll = TRUE, 
##     noise = TRUE, noise.factor = 1, force = FALSE, robMethod = "MM", 
##     force.mixed = TRUE, mi = 1, addMixedFactors = FALSE, trace = FALSE, 
##     init.method = "kNN", modelFormulas = NULL, multinom.method = "multinom", 
##     imp_var = TRUE, imp_suffix = "imp") 
## NULL

Rangkuman Algoritma EM adalah pendekatan komputasi untuk menemukan solusi untuk penaksir kemungkinan maksimum. Pada dasarnya, algoritma EM terdiri dari dua langkah, yaitu langkah E untuk estimasi parameter dan langkah M untuk maksimisasi sesuai dengan parameter yang sebenarnya. Algoritma ini biasanya konvergen dengan cepat dan diaplikasikan di banyak bidang.

Pada bab ini kita melihat penerapannya di dua area, yaitu pengelompokan dan imputasi nilai yang hilang. Pengelompokan adalah masalah NP-hard; secara umum, kita tidak dapat menemukan solusi bentuk tertutup yang tepat dalam waktu yang masuk akal. Oleh karena itu, algoritma EM diperlukan untuk menemukan solusi yang baik secara interaktif. Dalam pengelompokan, algoritma EM diimplementasikan untuk algoritma pengelompokan k-means, tetapi juga (tidak ditunjukkan dalam bab ini) untuk pengelompokan berbasis model dan untuk model campuran secara umum.

Nilai yang hilang sering terjadi dalam kumpulan data dalam praktiknya. Ilmuwan data mungkin adalah orang-orang yang pekerjaan utamanya adalah melakukan pra-pemrosesan data, sehingga mereka juga harus memperhitungkan nilai yang hilang. Kami melihat bahwa algoritma EM adalah alat utama untuk tugas ini.