mclust パッケージメモ

# install.packages('mclust')
library(mclust)
# load sample data
data(iris)
# fix column names
names(iris)<-tolower(names(iris))
# calc BIC
iris.clust <- mclustBIC(iris[,-5])
plot(iris.clust) # plot の結果からは2のmodelVEVがよさげ

plot of chunk mclust_data 

getClusteredData<-function(iris){
    # set model and data
    iris.hc        <- hc("VVV",iris[,-5]) # VEVモデルは非対応？VVVで。remove species column
    # num of clusters
    iris.cl        <- hclass(iris.hc,3)
    iris.result    <- cbind(iris,cl=iris.cl[,1])
    iris.result$cl <- factor(iris.result$cl)    
    return(iris.result)
}
iris.result <- getClusteredData(iris)
str(iris.result) 

## 'data.frame':    150 obs. of  6 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 ...
##  $ cl          : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...

# cross validation
ftable(iris.result[,c(5,6)])/50

##            cl    1    2    3
## species                     
## setosa        1.00 0.00 0.00
## versicolor    0.00 1.00 0.00
## virginica     0.00 0.28 0.72

library(GGally)
select.cols <- c("sepal.length","sepal.width","petal.length","petal.width")
# グラフのプロット　
ggpairs(
    title="実測グループ"
    ,data=iris
    ,columns=select.cols
    ,col="species")

plot of chunk mclust_plot 

ggpairs(
    title="予測グループ"
    ,data=iris.result
    ,columns=select.cols
    ,col="cl")

plot of chunk mclust_plot

参考リンク