R で可視化〜1つ1つの変数の分布の確認〜

2021年10月10日2022年1月12日R,Rで可視化,データサイエンス

ggplot2 を使って、1つ1つの変数の分布の確認を色々な方法でやってみます。

データは palmerpenguins を使います。このデータは Palmer Archipelago にある3つの島で観測された、3種類のペンギンの身長、体重などのデータです。

基本的な確認：データの形、型、欠損の有無

“palmerpenguins" を入手し、基本的な情報を見ていきます。また、特に言及しませんがパッケージは必要に応じてインストールをお願いします。

library(dplyr); library(tidyverse); library(ggplot2) # 毎回使うパッケージ
library(palmerpenguins)
print(penguins, n = 5, width = Inf) # 5行、全列表示

# 実行結果
# A tibble: 344 x 8
  species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g sex     year
  <fct>   <fct>              <dbl>         <dbl>             <int>       <int> <fct>  <int>
1 Adelie  Torgersen           39.1          18.7               181        3750 male    2007
2 Adelie  Torgersen           39.5          17.4               186        3800 female  2007
3 Adelie  Torgersen           40.3          18                 195        3250 female  2007
4 Adelie  Torgersen           NA            NA                  NA          NA NA      2007
5 Adelie  Torgersen           36.7          19.3               193        3450 female  2007
# … with 339 more rows

サンプルサイズが344、変数の数が8のクロスセクションデータであることがわかります。欠損値がどれほどあるのかを visdat::vis_miss() を使ってみてみます。

library(visdat)
vis_miss(x = penguins)

目視ですが、「欠損はペンギンの体型に関する測定値と性別で生じる」「体型に関する測定値はまとめて欠損する」「体型に関する測定値が欠損すれば性別も欠損する」といったことが伺えます。

本来であれば欠損値に対する対処はよく考えて行われるべきですが、今回はあくまで可視化が目的なので省略し、欠損を含むサンプルは削除します。

data = penguins %>% na.omit()

可視化

ggplot2 を使って、各変数がどんな分布をしているのかをみていきます。まず、ggplot() の初登場ということで簡単な図を作るために、分類対象となりそうな species と体に関する定量的な情報の一つである body_mass_g のそれぞれの分布を見てみます。独立に作成した一つ一つの図を一度に表示するには gridExtra パッケージが便利です。

library(gridExtra)
# 1つ目の図、ペンギンの種類の棒グラフ
fig1 = ggplot(data = data, aes(x = species)) + 
  geom_bar()
# 2つ目の図、体重のヒストグラム
fig2 = ggplot(data = data, aes(x = body_mass_g)) + 
  geom_histogram()

layout1 = matrix(1:2, nrow = 1)
fig_all = grid.arrange(fig1, fig2, layout_matrix = layout1)

ペンギンの種類を目的変数とするとするならば、極端な不均衡データとは言えないくらいの分布です。体重に関しては対数正規分布のような、歪度が正であるような分布であるように見えます。次に、body_mass_g を species で色分けして図示してみます。あるカテゴリによって色や点の形を変えたい場合は aes() を使って設定します。下図の場合、ビンごとの枠の色（col）と塗りつぶしの色（fill）を species ごとに設定しています。

ggplot(data = data, aes(x = body_mass_g)) + 
  geom_histogram(aes(col = species, fill = species), position = 'identity', alpha = 0.6)

Adelie ペンギンと Chinstrap ペンギンが同じような平均を持ち、Gentoo ペンギンだけ比較的重いことがわかります。全体としては対数正規分布のように見えましたが、種類ごとにみるとそれぞれ顕著な歪みはないことがわかります。もし species を分類する問題を考えるのであれば、少なくとも Gentoo かどうかに関しては体重が大きく貢献できそうです。このように、カテゴリ別に分布をみるだけで得られる情報量は大きく異なります。ちなみに、個別の頻度ではなく積み上げのヒストグラムとして表示することも可能です（position = 'stack'）。

同じような主旨の可視化を、ヒストグラム（相対頻度とカーネル密度推定による密度）、箱ひげ図、バイオリンプロット、散布図でもやってみます。なお、散布図のところで x = 1:nrow(data) とすることで x 軸の値を連番にし、データの順番通りにプロットしています。

# ヒストグラム（相対頻度とカーネル密度推定による密度）
fig_hist = ggplot(data = data, aes(x = body_mass_g, y = ..density..)) + 
  geom_histogram(aes(col = species, fill = species), position = 'identity', alpha = 0.6) +
  geom_density()
# 箱ひげ図
fig_box = ggplot(data = data, aes(x = species, y = body_mass_g)) + 
  geom_boxplot(aes(col = species))
# バイオリンプロット
fig_violin = ggplot(data = data, aes(x = species, y = body_mass_g)) + 
  geom_violin(aes(col = species))
# 散布図
fig_scatter = ggplot(data = data, aes(x = 1:nrow(data), y = body_mass_g)) + 
  geom_point(aes(col = species))

layout2 = matrix(1:4, nrow = 2)
fig_all = grid.arrange(fig_hist, fig_box, fig_violin, fig_scatter, layout_matrix = layout2)

分布を確認する手段にも様々あり、それぞれに長短あります。例えば箱ひげ図では四分位点が分かりますがどんな形をしているかが見えません、バイオリンプロットでは滑らかにされた相対頻度として分布の形が見えますが、外れ値や分布の詳細な凸凹具合はほとんどなくなります。また、できれば別々に図を表示して考えるよりも、より少ない図でデータを理解したいものです。

これに対する策として、異なる図を重ね書きするということをやってみます。今回はバイオリンプロット、箱ひげ図、サンプルの点列（ランダムにぶれさせたもの）を重ね書きしてみます。

ggplot(data = data, aes(x = species, y = body_mass_g)) + 
  geom_violin(aes(col = species, fill = species), alpha = 0.5) +
  geom_boxplot(width = 0.1, alpha = 0.5) + 
  geom_jitter(size = 0.3, width = 0.1, alpha = 0.4)

重ね書きといっても煩雑にならず、スッキリまとまった図になったのではないでしょうか。ちなみに「サンプルの点列（ランダムにぶれさせたもの）」は geom_jitter() で書かれていますが、これは geom_point(position = 'jitter') と同じで、サンプルがどの値のところで集中しているかをみやすくするために左右に適当にぶれさせて表示させるものです。なので、左右の位置に意味はありません。

可能であればこのようなグラフを全ての変数に対して作成し、まとめて表示させて見比べたいものです。次はこの「全ての変数の分布をまとめて表示」をやってみみたいと思います。

まずはまとめて図示するためにデータを加工する必要があります。ざっくりいうと、「プロットする値（value）」「変数名（variable）」「カテゴリ（id）」という3つの列に分けてデータを整理する必要があります。実際にデータを見た方が理解が早いので、reshape2::melt() を使ってデータを加工するところと結果を見てみます。なお、データは分類基準にするカテゴリ（今回は species）と連続値に絞っています。

library(reshape2)
data_gg <- melt(data[,c(1,3:6)], id = 'species') # species と連続値に限定
head(data_gg); tail(data_gg)

# 実行結果
 species       variable value
1  Adelie bill_length_mm  39.1
2  Adelie bill_length_mm  39.5
3  Adelie bill_length_mm  40.3
4  Adelie bill_length_mm  36.7
5  Adelie bill_length_mm  39.3
6  Adelie bill_length_mm  38.9
       species    variable value
1327 Chinstrap body_mass_g  3650
1328 Chinstrap body_mass_g  4000
1329 Chinstrap body_mass_g  3400
1330 Chinstrap body_mass_g  3775
1331 Chinstrap body_mass_g  4100
1332 Chinstrap body_mass_g  3775

このように値、変数名、カテゴリに分けた縦長のデータを作ることで「species（id）でグルーピングした値（value）を変数（variable）ごとに図示」ということがやりやすくなります。

データが整ったら可視化の準備完了です。species ごとに連続値に該当する全ての変数を可視化してみます。ある可視化を変数ごとに独立に行って同時に表示する際は facet_wrap() を使います。~variable は「variable ごとに分割する」という指示で、scales = 'free' とすることで変数ごとに y 軸のスケールを可変にしています。逆に、全ての図に対して一定のスケールを使用する場合は scales = 'fixed' と指定します。

ggplot(data = data_gg, aes(x = species, y = value)) +
  facet_wrap(~variable, scales = 'free') +
  geom_violin(aes(col = species, fill = species), alpha = 0.5) +
  geom_boxplot(width = 0.1, alpha = 0.5) + 
  geom_jitter(size = 0.2, width = 0.2, alpha = 0.4)