R {ggplot2} の散布図に凸包 / 確率楕円を描きたい

小ネタ。{ggplot2} でグループ別の散布図を描くときに、ちょっと飾り付けをしてグループをわかりやすくしたい。

凸包 (Convex)

最初にベースとなる散布図を描く。

library(dplyr)
library(ggplot2)

df <- iris

p <- ggplot(df, aes(x = Petal.Width, y = Petal.Length)) + 
  geom_point()
p

まずは 散布図全体について凸包をとる。ある点の集合の凸包は、 grDevices::chull で計算できる。chull は凸な点の index を返すので、この返り値に含まれるデータのみをフィルタして geom_polygon に渡せばよい。

chull(df[c('Petal.Width', 'Petal.Length')])
#  [1]  44  17  23  14  33  25 135 123 119 110 145 115

hulls <- df[chull(df[c('Petal.Width', 'Petal.Length')]), ]
hulls
#     Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
# 44           5.0         3.5          1.6         0.6    setosa
# 17           5.4         3.9          1.3         0.4    setosa
# 23           4.6         3.6          1.0         0.2    setosa
# 14           4.3         3.0          1.1         0.1    setosa
# 33           5.2         4.1          1.5         0.1    setosa
# 25           4.8         3.4          1.9         0.2    setosa
# 135          6.1         2.6          5.6         1.4 virginica
# 123          7.7         2.8          6.7         2.0 virginica
# 119          7.7         2.6          6.9         2.3 virginica
# 110          7.2         3.6          6.1         2.5 virginica
# 145          6.7         3.3          5.7         2.5 virginica
# 115          5.8         2.8          5.1         2.4 virginica

p + geom_polygon(data = hulls, alpha = 0.2)

凸包をグループ別に描画したい場合は、{dplyr} で各グループへ chull を適用する。

p <- ggplot(df, aes(x = Petal.Width, y = Petal.Length, colour = Species)) + 
  geom_point()

hulls <- df %>%
  dplyr::group_by(Species) %>%
  dplyr::do(.[chull(.[c('Petal.Width', 'Petal.Length')]), ])
hulls
# Source: local data frame [22 x 5]
# Groups: Species
# 
#    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
# 1           5.4         3.9          1.3         0.4     setosa
# 2           4.6         3.6          1.0         0.2     setosa
# 3           4.3         3.0          1.1         0.1     setosa
# 4           5.2         4.1          1.5         0.1     setosa
# 5           4.8         3.4          1.9         0.2     setosa
# 6           5.1         3.8          1.9         0.4     setosa
# 7           5.0         3.5          1.6         0.6     setosa
# 8           5.1         2.5          3.0         1.1 versicolor
# 9           4.9         2.4          3.3         1.0 versicolor
# 10          5.8         2.7          4.1         1.0 versicolor
# 11          6.1         2.8          4.7         1.2 versicolor
# 12          6.0         2.7          5.1         1.6 versicolor
# 13          6.7         3.0          5.0         1.7 versicolor
# 14          5.9         3.2          4.8         1.8 versicolor
# 15          5.8         2.8          5.1         2.4  virginica
# 16          4.9         2.5          4.5         1.7  virginica
# 17          6.0         2.2          5.0         1.5  virginica
# 18          6.1         2.6          5.6         1.4  virginica
# 19          7.7         2.8          6.7         2.0  virginica
# 20          7.7         2.6          6.9         2.3  virginica
# 21          7.2         3.6          6.1         2.5  virginica
# 22          6.7         3.3          5.7         2.5  virginica

p + geom_polygon(data = hulls, alpha = 0.2)

参考 以下リンクにある plyr::ddply を使うバージョンを参考にした。

• How to draw neat polygons around scatterplot regions in ggplot2 - Cross Validated

確率楕円 (Probability Ellipse)

こっちは ggplot2 1.0.0 以降なら stat_ellipse 一発なので簡単。

p + stat_ellipse()

クラスタリング結果への凸包 / 確率楕円の描画

もともとやりたかったのは {ggplot2} + {ggfortify} でクラスタリング結果を autoplot するときに凸包 / 確率楕円を描画すること。これを autoplot 時の frame オプションで指定できるようにした。

ついでに {cluster} パッケージの非階層クラスタリング法 cluster::clara, cluster::fanny, cluster::pam も autoplot できるようにした。

# library(devtools)
# install_github('sinhrks/ggfortify')

library(ggplot2)
library(ggfortify)
df <- iris[-5]

autoplot(kmeans(df, 3), original = iris, frame = TRUE)

library(cluster)
autoplot(clara(df, 3), frame = TRUE, frame.type = 'norm')

ほか、細かい使い方はこちら。

• RPubs - Plotting PCA/clustering results using ggplot2 and ggfortify

クラスタリング以外の autoplot についてはこちら。

• 線形回帰 / 一般化線形モデル
• 時系列
• 時系列関連の統計
• 多次元尺度構成法
• 生存曲線
• 確率分布

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集