Analisis cluster di R: menentukan jumlah cluster yang optimal

428

Menjadi pemula di R, saya tidak begitu yakin bagaimana memilih jumlah cluster terbaik untuk melakukan analisis k-means. Setelah merencanakan subkumpulan data di bawah ini, berapa banyak kluster yang akan sesuai? Bagaimana saya bisa melakukan analisis dendro klaster?

n = 1000
kk = 10    
x1 = runif(kk)
y1 = runif(kk)
z1 = runif(kk)    
x4 = sample(x1,length(x1))
y4 = sample(y1,length(y1)) 
randObs <- function()
{
  ix = sample( 1:length(x4), 1 )
  iy = sample( 1:length(y4), 1 )
  rx = rnorm( 1, x4[ix], runif(1)/8 )
  ry = rnorm( 1, y4[ix], runif(1)/8 )
  return( c(rx,ry) )
}  
x = c()
y = c()
for ( k in 1:n )
{
  rPair  =  randObs()
  x  =  c( x, rPair[1] )
  y  =  c( y, rPair[2] )
}
z <- rnorm(n)
d <- data.frame( x, y, z )
pengguna2153893
sumber
4
Jika Anda tidak sepenuhnya terhubung dengan kmeans, Anda bisa mencoba algoritma pengelompokan DBSCAN, tersedia dalam fpcpaket. Memang benar, Anda kemudian harus menetapkan dua parameter ... tetapi saya telah menemukan bahwa fpc::dbscankemudian melakukan pekerjaan yang cukup baik secara otomatis menentukan jumlah cluster yang baik. Plus itu sebenarnya dapat menghasilkan satu cluster jika data itu memberi tahu Anda - beberapa metode dalam jawaban luar biasa @ Ben tidak akan membantu Anda menentukan apakah k = 1 sebenarnya yang terbaik.
Stephan Kolassa
Lihat juga stats.stackexchange.com/q/11691/478
Richie Cotton

Jawaban:

1020

Jika pertanyaan Anda adalah how can I determine how many clusters are appropriate for a kmeans analysis of my data?, maka berikut adalah beberapa opsi. The artikel wikipedia pada penentuan jumlah cluster memiliki review yang baik tentang beberapa metode ini.

Pertama, beberapa data yang dapat direproduksi (data dalam Q itu ... tidak jelas bagi saya):

n = 100
g = 6 
set.seed(g)
d <- data.frame(x = unlist(lapply(1:g, function(i) rnorm(n/g, runif(1)*i^2))), 
                y = unlist(lapply(1:g, function(i) rnorm(n/g, runif(1)*i^2))))
plot(d)

masukkan deskripsi gambar di sini

Satu . Cari tikungan atau siku dalam jumlah plot scree jumlah kuadrat (SSE). Lihat http://www.statmethods.net/advstats/cluster.html & http://www.mattpeeples.net/kmeans.html untuk informasi lebih lanjut. Lokasi siku di plot yang dihasilkan menunjukkan jumlah cluster yang cocok untuk kmeans:

mydata <- d
wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))
  for (i in 2:15) wss[i] <- sum(kmeans(mydata,
                                       centers=i)$withinss)
plot(1:15, wss, type="b", xlab="Number of Clusters",
     ylab="Within groups sum of squares")

Kami dapat menyimpulkan bahwa 4 cluster akan ditunjukkan dengan metode ini: masukkan deskripsi gambar di sini

Dua . Anda dapat melakukan partisi di sekitar medoid untuk memperkirakan jumlah cluster menggunakan pamkfungsi dalam paket fpc.

library(fpc)
pamk.best <- pamk(d)
cat("number of clusters estimated by optimum average silhouette width:", pamk.best$nc, "\n")
plot(pam(d, pamk.best$nc))

masukkan deskripsi gambar di sini masukkan deskripsi gambar di sini

# we could also do:
library(fpc)
asw <- numeric(20)
for (k in 2:20)
  asw[[k]] <- pam(d, k) $ silinfo $ avg.width
k.best <- which.max(asw)
cat("silhouette-optimal number of clusters:", k.best, "\n")
# still 4

Tiga . Kriteria Calinsky: Pendekatan lain untuk mendiagnosis berapa banyak cluster yang sesuai dengan data. Dalam hal ini kami mencoba 1 hingga 10 grup.

require(vegan)
fit <- cascadeKM(scale(d, center = TRUE,  scale = TRUE), 1, 10, iter = 1000)
plot(fit, sortg = TRUE, grpmts.plot = TRUE)
calinski.best <- as.numeric(which.max(fit$results[2,]))
cat("Calinski criterion optimal number of clusters:", calinski.best, "\n")
# 5 clusters!

masukkan deskripsi gambar di sini

Empat . Tentukan model yang optimal dan jumlah cluster berdasarkan Kriteria Informasi Bayesian untuk maksimalisasi-harapan, diinisialisasi oleh pengelompokan hierarkis untuk model campuran Gaussian berparameterisasi

# See http://www.jstatsoft.org/v18/i06/paper
# http://www.stat.washington.edu/research/reports/2006/tr504.pdf
#
library(mclust)
# Run the function to see how many clusters
# it finds to be optimal, set it to search for
# at least 1 model and up 20.
d_clust <- Mclust(as.matrix(d), G=1:20)
m.best <- dim(d_clust$z)[2]
cat("model-based optimal number of clusters:", m.best, "\n")
# 4 clusters
plot(d_clust)

masukkan deskripsi gambar di sini masukkan deskripsi gambar di sini masukkan deskripsi gambar di sini

Lima . Pengelompokan afinitas (AP), lihat http://dx.doi.org/10.1126/science.1136800

library(apcluster)
d.apclus <- apcluster(negDistMat(r=2), d)
cat("affinity propogation optimal number of clusters:", length(d.apclus@clusters), "\n")
# 4
heatmap(d.apclus)
plot(d.apclus, d)

masukkan deskripsi gambar di sini masukkan deskripsi gambar di sini

Enam . Statistik Gap untuk Memperkirakan Jumlah Cluster. Lihat juga beberapa kode untuk output grafis yang bagus . Mencoba 2-10 cluster di sini:

library(cluster)
clusGap(d, kmeans, 10, B = 100, verbose = interactive())

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 100)  [one "." per sample]:
.................................................. 50 
.................................................. 100 
Clustering Gap statistic ["clusGap"].
B=100 simulated reference sets, k = 1..10
 --> Number of clusters (method 'firstSEmax', SE.factor=1): 4
          logW   E.logW        gap     SE.sim
 [1,] 5.991701 5.970454 -0.0212471 0.04388506
 [2,] 5.152666 5.367256  0.2145907 0.04057451
 [3,] 4.557779 5.069601  0.5118225 0.03215540
 [4,] 3.928959 4.880453  0.9514943 0.04630399
 [5,] 3.789319 4.766903  0.9775842 0.04826191
 [6,] 3.747539 4.670100  0.9225607 0.03898850
 [7,] 3.582373 4.590136  1.0077628 0.04892236
 [8,] 3.528791 4.509247  0.9804556 0.04701930
 [9,] 3.442481 4.433200  0.9907197 0.04935647
[10,] 3.445291 4.369232  0.9239414 0.05055486

Inilah output dari implementasi statistik gap Edwin Chen: masukkan deskripsi gambar di sini

Tujuh . Anda juga mungkin merasa berguna untuk mengeksplorasi data Anda dengan clustergram untuk memvisualisasikan penetapan cluster, lihat http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r- kode / untuk lebih jelasnya.

Delapan . The NbClust paket menyediakan 30 indeks untuk menentukan jumlah cluster dalam kumpulan data.

library(NbClust)
nb <- NbClust(d, diss=NULL, distance = "euclidean",
        method = "kmeans", min.nc=2, max.nc=15, 
        index = "alllong", alphaBeale = 0.1)
hist(nb$Best.nc[1,], breaks = max(na.omit(nb$Best.nc[1,])))
# Looks like 3 is the most frequently determined number of clusters
# and curiously, four clusters is not in the output at all!

masukkan deskripsi gambar di sini

Jika pertanyaan Anda adalah how can I produce a dendrogram to visualize the results of my cluster analysis, maka Anda harus mulai dengan ini: http://www.statmethods.net/advstats/cluster.html http://www.r-tutor.com/gpu-computing/clustering/hierarchical-cluster-analysis http://gastonsanchez.wordpress.com/2012/10/03/7-ways-to-plot-dendrograms-in-r/ Dan lihat di sini untuk metode yang lebih eksotis: http://cran.r-project.org/ web / views / Cluster.html

Berikut ini beberapa contoh:

d_dist <- dist(as.matrix(d))   # find distance matrix 
plot(hclust(d_dist))           # apply hirarchical clustering and plot

masukkan deskripsi gambar di sini

# a Bayesian clustering method, good for high-dimension data, more details:
# http://vahid.probstat.ca/paper/2012-bclust.pdf
install.packages("bclust")
library(bclust)
x <- as.matrix(d)
d.bclus <- bclust(x, transformed.par = c(0, -50, log(16), 0, 0, 0))
viplot(imp(d.bclus)$var); plot(d.bclus); ditplot(d.bclus)
dptplot(d.bclus, scale = 20, horizbar.plot = TRUE,varimp = imp(d.bclus)$var, horizbar.distance = 0, dendrogram.lwd = 2)
# I just include the dendrogram here

masukkan deskripsi gambar di sini

Juga untuk data dimensi tinggi adalah pvclustpustaka yang menghitung nilai-p untuk pengelompokan hierarkis melalui resampling bootstrap multiskala. Inilah contoh dari dokumentasi (tidak akan bekerja pada data dimensi rendah seperti pada contoh saya):

library(pvclust)
library(MASS)
data(Boston)
boston.pv <- pvclust(Boston)
plot(boston.pv)

masukkan deskripsi gambar di sini

Apakah ada yang membantu?

Ben
sumber
Untuk dendogram terakhir (Cluster Dendogram dengan AU / BP) kadang-kadang lebih mudah untuk menggambar persegi panjang di sekitar kelompok dengan nilai-p yang relatif tinggi: pvrect (fit, alpha = 0,95)
Igor Elbert
Inilah yang saya cari. Saya baru mengenal R dan perlu waktu lama untuk menemukan ini. Terima kasih @Ben untuk menjawab dengan detail seperti itu. Bisakah Anda membimbing saya ke mana saya dapat menemukan logika di balik masing-masing metode ini, seperti apa metrik atau kriteria yang mereka gunakan untuk menentukan jumlah cluster yang optimal, atau bagaimana masing-masing dari mereka berbeda satu sama lain. Bos saya ingin saya mengatakannya, sehingga kami dapat memutuskan metode mana yang akan digunakan. Terima kasih sebelumnya.
nasia jaffri
1
@Alexandr Blekh Anda juga dapat mencoba mengubah metode grafis apa pun menjadi analitis. Misalnya, saya menggunakan metode "siku" (pertama disebutkan dalam jawaban), tetapi cobalah untuk menemukannya secara analitis. Titik siku bisa menjadi titik dengan kelengkungan maksimum. Untuk data diskrit, ini adalah titik dengan perbedaan pusat urutan kedua maksimum (turunan urutan analog ke maks. Kedua untuk data kontinu). Lihat stackoverflow.com/a/4473065/1075993 dan stackoverflow.com/q/2018178/1075993 . Saya rasa metode grafis lain juga dapat dikonversi menjadi analitis.
Andrey Sapegin
1
@AndreySapegin: Saya bisa, tetapi: 1) terus terang, saya tidak menganggapnya sebagai solusi yang elegan (IMHO, dalam banyak kasus, metode visual harus tetap visual, sedangkan yang analitis harus tetap analitis); 2) Saya telah menemukan solusi analitik untuk ini, menggunakan satu atau beberapa Rpaket (ada di GitHub saya - Anda dapat melihatnya); 3) solusi saya tampaknya bekerja cukup baik, ditambah, sudah lama dan saya sudah menyelesaikan perangkat lunak disertasi saya, laporan disertasi (tesis) dan saat ini saya sedang mempersiapkan pertahanan :-). Bagaimanapun, saya sangat menghargai komentar dan tautan Anda. Semua yang terbaik!
Aleksandr Blekh
1
2,2 juta baris ada di dataset clustering saya saat ini. Tidak ada paket R ini yang berfungsi, saya harapkan. Mereka hanya mem-pop komputer saya dan kemudian jatuh dari pengalaman saya. Namun, sepertinya penulis tahu barang-barangnya untuk data kecil dan untuk kasus umum tanpa memperhatikan kapasitas perangkat lunak. Tidak ada poin yang dikurangi karena kerja bagus yang jelas dari penulis. Kalian hanya perlu tahu bahwa R lama biasa mengerikan pada 2,2 juta baris - coba sendiri jika Anda tidak percaya padaku. H2O membantu tetapi terbatas pada taman kecil berdinding bahagia.
Geoffrey Anderson
21

Sulit untuk menambahkan sesuatu jawaban yang terlalu rumit. Meskipun saya merasa kita harus menyebutkannya di identifysini, terutama karena @Ben menunjukkan banyak contoh dendrogram.

d_dist <- dist(as.matrix(d))   # find distance matrix 
plot(hclust(d_dist)) 
clusters <- identify(hclust(d_dist))

identifymemungkinkan Anda secara interaktif memilih cluster dari dendrogram dan menyimpan pilihan Anda ke daftar. Tekan Esc untuk meninggalkan mode interaktif dan kembali ke konsol R. Perhatikan, bahwa daftar itu berisi indeks, bukan rownames (berlawanan dengan cutree).

Matt Bannert
sumber
10

Untuk menentukan k-cluster optimal dalam metode clustering. Saya biasanya menggunakan Elbowmetode yang disertai oleh pemrosesan Paralel untuk menghindari pemulaan waktu. Kode ini dapat dicontoh seperti ini:

Metode siku

elbow.k <- function(mydata){
dist.obj <- dist(mydata)
hclust.obj <- hclust(dist.obj)
css.obj <- css.hclust(dist.obj,hclust.obj)
elbow.obj <- elbow.batch(css.obj)
k <- elbow.obj$k
return(k)
}

Menjalankan Elbow parallel

no_cores <- detectCores()
    cl<-makeCluster(no_cores)
    clusterEvalQ(cl, library(GMD))
    clusterExport(cl, list("data.clustering", "data.convert", "elbow.k", "clustering.kmeans"))
 start.time <- Sys.time()
 elbow.k.handle(data.clustering))
 k.clusters <- parSapply(cl, 1, function(x) elbow.k(data.clustering))
    end.time <- Sys.time()
    cat('Time to find k using Elbow method is',(end.time - start.time),'seconds with k value:', k.clusters)

Ini bekerja dengan baik.

VanThaoNguyen
sumber
2
Fungsi siku dan css berasal dari paket GMD: cran.r-project.org/web/packages/GMD/GMD.pdf
Rohan
6

Jawaban indah dari Ben. Namun saya terkejut bahwa metode Affinity Propagation (AP) telah di sini disarankan hanya untuk menemukan jumlah cluster untuk metode k-means, di mana secara umum AP melakukan pekerjaan yang lebih baik mengelompokkan data. Silakan lihat makalah ilmiah yang mendukung metode ini dalam Sains di sini:

Frey, Brendan J., dan Delbert Dueck. "Clustering dengan mengirimkan pesan antar titik data." sains 315.5814 (2007): 972-976.

Jadi jika Anda tidak bias terhadap k-means saya sarankan untuk menggunakan AP secara langsung, yang akan mengelompokkan data tanpa perlu mengetahui jumlah cluster:

library(apcluster)
apclus = apcluster(negDistMat(r=2), data)
show(apclus)

Jika jarak euclidean negatif tidak sesuai, maka Anda dapat menggunakan langkah-langkah kesamaan lainnya yang disediakan dalam paket yang sama. Misalnya, untuk persamaan berdasarkan korelasi Spearman, inilah yang Anda butuhkan:

sim = corSimMat(data, method="spearman")
apclus = apcluster(s=sim)

Harap dicatat bahwa fungsi-fungsi untuk kesamaan dalam paket AP hanya disediakan untuk kesederhanaan. Bahkan, fungsi apcluster () dalam R akan menerima matriks korelasi apa pun. Hal yang sama sebelumnya dengan corSimMat () dapat dilakukan dengan ini:

sim = cor(data, method="spearman")

atau

sim = cor(t(data), method="spearman")

tergantung pada apa yang ingin Anda klaster pada matriks Anda (baris atau cols).

zsram
sumber
6

Metode ini bagus, tetapi ketika mencoba mencari k untuk set data yang jauh lebih besar, ini bisa sangat lambat di R.

Solusi bagus yang saya temukan adalah paket "RWeka", yang memiliki implementasi efisien dari algoritma X-Means - versi tambahan dari K-Means yang memiliki skala yang lebih baik dan akan menentukan jumlah cluster yang optimal untuk Anda.

Pertama, Anda ingin memastikan bahwa Weka diinstal pada sistem Anda dan memiliki XMeans diinstal melalui alat manajer paket Weka.

library(RWeka)

# Print a list of available options for the X-Means algorithm
WOW("XMeans")

# Create a Weka_control object which will specify our parameters
weka_ctrl <- Weka_control(
    I = 1000,                          # max no. of overall iterations
    M = 1000,                          # max no. of iterations in the kMeans loop
    L = 20,                            # min no. of clusters
    H = 150,                           # max no. of clusters
    D = "weka.core.EuclideanDistance", # distance metric Euclidean
    C = 0.4,                           # cutoff factor ???
    S = 12                             # random number seed (for reproducibility)
)

# Run the algorithm on your data, d
x_means <- XMeans(d, control = weka_ctrl)

# Assign cluster IDs to original data set
d$xmeans.cluster <- x_means$class_ids
RDRR
sumber
6

Solusi sederhana adalah perpustakaan factoextra. Anda dapat mengubah metode pengelompokan dan metode untuk menghitung jumlah grup terbaik. Misalnya jika Anda ingin mengetahui jumlah cluster terbaik untuk k- berarti:

Data: mtcars

library(factoextra)   
fviz_nbclust(mtcars, kmeans, method = "wss") +
      geom_vline(xintercept = 3, linetype = 2)+
      labs(subtitle = "Elbow method")

Akhirnya, kita mendapatkan grafik seperti:

masukkan deskripsi gambar di sini

Cro-Magnon
sumber
2

Jawabannya bagus. Jika Anda ingin memberikan kesempatan kepada metode pengelompokan lain, Anda dapat menggunakan pengelompokan hierarkis dan melihat bagaimana data terpecah.

> set.seed(2)
> x=matrix(rnorm(50*2), ncol=2)
> hc.complete = hclust(dist(x), method="complete")
> plot(hc.complete)

masukkan deskripsi gambar di sini

Bergantung pada berapa kelas yang Anda butuhkan, Anda dapat memotong dendrogram sebagai;

> cutree(hc.complete,k = 2)
 [1] 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 2 1 1 1
[26] 2 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 1 1 1 1 1 1 1 2

Jika Anda mengetik, ?cutreeAnda akan melihat definisi. Jika kumpulan data Anda memiliki tiga kelas, itu akan menjadi sederhana cutree(hc.complete, k = 3). Setara dengan cutree(hc.complete,k = 2)adalah cutree(hc.complete,h = 4.9).

boyaronur
sumber
Saya lebih suka Bangsal lebih lengkap.
Chris