1 本講義の目的

  • 要約統計量の算出方法を学ぶ。
  • ggplotやplotlyなどを用いたデータの可視化の基礎を学ぶ。

2 パッケージの読み込み

pacmanを使うと、install.packages()やlibrary()を使わずにパッケージをまとめて読み込むことができる。

まずpackmanをインストールする。

install.packages("pacman")

今回は、readxlとtidyverseとplotlyとskimrを使う。以下のように書けば、まとめて読み込むことができる。

pacman::p_load(readxl, tidyverse, plotly, skimr)

3 ディレクトリの設定

（RStudio CloudではなくローカルPCで作業する場合）作業ディレクトリを設定する。 プロジェクト機能を活用できるのであれば、そちらを活用する。

setwd(自分のディレクトリ) # ディレクトリの設定
getwd() # ディレクトリの確認

4 例示用データの用意と前処理

前回同様、今回も test_scores.xlsx のデータを使う

df <- read_excel("test_scores.xlsx") # エクセル読み込み
head(df)
## # A tibble: 6 × 5
##   クラス 名前   数学  英語  国語
##   <chr>  <chr> <dbl> <dbl> <dbl>
## 1 C      安藤     67    78    45
## 2 A      藤村     55    73    72
## 3 B      大泉     47    55    73
## 4 B      木村     78    71    85
## 5 A      鈴井     51    74    61
## 6 C      望月     79    84    75

今回は、dplyr::rename()によって変数名を英語に変換してから使う。

df <- df %>% dplyr::rename(class = クラス,
                           name = 名前,
                           math = 数学,
                           english = 英語,
                           japanese = 国語)
head(df)
## # A tibble: 6 × 5
##   class name   math english japanese
##   <chr> <chr> <dbl>   <dbl>    <dbl>
## 1 C     安藤     67      78       45
## 2 A     藤村     55      73       72
## 3 B     大泉     47      55       73
## 4 B     木村     78      71       85
## 5 A     鈴井     51      74       61
## 6 C     望月     79      84       75

5 データの分布を要約する

5.1 平均・分散など

データの分布を要約するための指標（要約統計量）の代表的なものが平均（mean）と分散（variance）である。

  • 平均（mean）：\(n\)個のデータ\(x_i\)の総和を、サンプルサイズ\(n\)で割ったもの。

\[ \bar{x} = \frac{1}{n}\sum_{i=1}^n x_i \]

  • 分散（variance）：平均からの差の2乗和の平均

\[ \sigma^2=\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2 \]

ただし、Rのvar()関数の算出する値は、母集団の分散（母分散）の不偏推定量としての分散すなわち不偏分散

\[ s^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2 \]

である。

# 不偏分散
var(df$english)
## [1] 333.1481

分散の単位は元の単位の2乗であるため、それを元の単位に戻すために平方根をとったものが標準偏差（standard deviation）である。

# 不偏標準偏差
sd(df$english)
## [1] 18.25234

5.2 要約統計量

5.2.1 summary()

summary()関数は、データフレームの各行に関する代表値などの情報を取得することができる関数である。

summary(df)
##     class               name                math           english     
##  Length:40          Length:40          Min.   : 42.00   Min.   :42.00  
##  Class :character   Class :character   1st Qu.: 51.75   1st Qu.:50.75  
##  Mode  :character   Mode  :character   Median : 65.00   Median :65.50  
##                                        Mean   : 66.58   Mean   :68.33  
##                                        3rd Qu.: 80.50   3rd Qu.:82.25  
##                                        Max.   :100.00   Max.   :99.00  
##     japanese     
##  Min.   : 40.00  
##  1st Qu.: 51.00  
##  Median : 67.50  
##  Mean   : 66.88  
##  3rd Qu.: 79.50  
##  Max.   :100.00

classとnameは文字列（character）型なので、Length（行数、標本数）が40であることくらいしか表示されない。

english、japanese、mathは整数（integer）型なので、代表値などの要約統計量（summary statistics）が表示されている。

それぞれの統計量の意味は

  • Min.が最小値
  • 1st Qu.が第一四分位数（1st quantile）：「データを小さいものから順に並べた時に標本数の4分の1の位置（下位25%）にくる値」
  • Medianが中央値（median）：「データを小さい順に並べた時に標本数の2分の1の位置にくる値」
  • Meanが平均値（mean）：「データの総計をデータの個数で割ったもの」\(\sum_{i=1}^n x_i/n\)
  • 3rd Qu.が第三四分位数（3rd quantile）：「データを小さい順に並べた時に標本数の4分の3の位置（75%）にくる値」
  • Max.が最大値

である。これらのうちMean以外の５つの統計量をまとめて五数要約という。

平均値は外れ値（極端に高かったり低かったりする値）の影響を受けやすいため、これらの統計量を見ることで、中央値と平均値が大きく異なるデータは外れ値があったりデータの分布に歪みがあることが推測できるわけである。

5.2.2 skimr::skim()

skimrのskim()を使うと、小さなヒストグラムつきの要約統計量の表を出力できる。欠損値の数なども表示してくれる。

skimr::skim(df)
Data summary
Name df
Number of rows 40
Number of columns 5
_______________________
Column type frequency:
character 2
numeric 3
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
class 0 1 1 1 0 3 0
name 0 1 1 3 0 38 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
math 0 1 66.58 18.13 42 51.75 65.0 80.50 100 ▇▆▅▃▅
english 0 1 68.33 18.25 42 50.75 65.5 82.25 99 ▇▅▂▆▅
japanese 0 1 66.88 17.17 40 51.00 67.5 79.50 100 ▇▃▇▅▃

6 2変数の関係を要約する

6.1 相関係数

相関係数（correlation coefficient）は２つの量的変数の直線的な関係の強さを測るもので、次の式で定義される

\[ \begin{align} x\text{と}y\text{の相関係数}r &= \frac{x\text{と}y\text{の共分散}}{x\text{の標準偏差}\times y\text{の標準偏差}} \\ &= \frac{\frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2} \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i - \bar{y})^2}} \end{align} \]

Rではcor()で計算できる。相関係数は-1から+1の範囲の値をとり、絶対値で1に近いほど関係が強いことを示す。

cor(x = df$english, y = df$math)
## [1] 0.160109

6.2 相関行列

相関係数が入った行列を相関行列という

cor(df[,3:5])
##              math    english   japanese
## math     1.000000 0.16010896 0.22593799
## english  0.160109 1.00000000 0.04096556
## japanese 0.225938 0.04096556 1.00000000

7 描画用パッケージ

グラフの描画で役立つパッケージを紹介する。

7.1 ggplot2

{ggplot2}というグラフの描画用パッケージは綺麗なグラフを描くことができる。

ggplotは次のように２つ以上の手順をとり、２つ以上の関数を+でつなげていき、レイヤーを重ねるように描画していく。

  1. キャンバスを用意する（ggplot(データフレーム, aes(x, y))）
  2. グラフを用意する（geom_histogram()などgeom_なんたらという関数）
  3. グラフの装飾のための関数を使う場合は、さらにレイヤーを重ねていく
# 棒グラフをつくり、gに代入する
g <- ggplot(df, aes(x = `class`)) + # dfのデータでキャンバスを作る（`で囲っているのは関数のclass()と名前が被っているため）
  geom_bar() # キャンバスに棒グラフgeom_bar()を重ねる

g # 表示

（参考） ggplot2のチートシートは、公式サイトの「Japanese Translations – 日本語翻訳」で見ることができる。

7.2 plotly

{plotly}パッケージを使うと、簡単なインタラクティブなグラフ（カーソルを合わせると情報が表示されたり、表示情報を変更したりできるグラフ）を描くことができる。

ggplotly()関数を使うと、{ggplot2}で描いたグラフをインタラクティブなものに変換できる。

plotly::ggplotly(g)   # ggplot2のグラフをplotlyのグラフへ変換

8 1変数のグラフ

8.1 ヒストグラム

ヒストグラムは量的変数の分布を確認するには極めて有用なグラフである。

# ヒストグラム
ggplot(df, aes(x = english)) +
  geom_histogram(bins = 10) # binsは棒の数。指定しない場合はbins=30になる

なおR本体の機能でヒストグラムを描く場合はhist()を使う

# ヒストグラム
hist(x = df[["math"]])

8.2 棒グラフ

# ggplotの棒グラフ
g <- ggplot(df, aes(x = `class`)) +
  geom_bar()

#グラフ表示
g

#plolyで表示
plotly::ggplotly(g)

9 2変数のグラフ

9.1 散布図

2変数でグラフを描く場合、キャンバスにx（グラフ横軸）とy（グラフ縦軸）の両方を指定する。

そして散布図を描く場合はgeom_point()を+でつなげればよい

# ggplot
g <- ggplot(df, aes(x = english, y = japanese)) +
  geom_point() +
  xlab("English") + # X軸のラベルを変更
  ylab("Japanese") # y軸のラベルを変更

# 散布図表示
g

# plotlyで表示
plotly::ggplotly(g)

回帰直線を加えるには、stat_smooth()を使う。

ggplot(df, aes(x = english, y = japanese)) +
  geom_point() +
  stat_smooth(method = "lm", se =FALSE, colour = "blue") # lm: 最小二乗法、信頼区間は表示させない(FALSE)、色は青(blue)

9.2 折れ線グラフ

折れ線グラフはgeom_line()で描くことができる。

Excelで言うところの「マーカー付き折れ線グラフ」にしたい場合はgeom_point()を重ねればよい。

ggplot(df, aes(x = english, y = japanese)) +
  geom_line() + # 折れ線
  geom_point()  # マーカー

（参考）アメリカのGNP

上では説明の簡単のため散布図と同じデータを使ったが、折れ線グラフは本来は時系列データに使うのが望ましい。

そこで、以下ではRにあらかじめ収録されている練習用のデータセットのひとつであるlongleyを使った例をのせる。

このデータは1947~62年のアメリカのGNPや雇用者数などが収録されている（詳細はhelpを参照してほしい）。

head(longley)
##      GNP.deflator     GNP Unemployed Armed.Forces Population Year Employed
## 1947         83.0 234.289      235.6        159.0    107.608 1947   60.323
## 1948         88.5 259.426      232.5        145.6    108.632 1948   61.122
## 1949         88.2 258.054      368.2        161.6    109.773 1949   60.171
## 1950         89.5 284.599      335.1        165.0    110.929 1950   61.187
## 1951         96.2 328.975      209.9        309.9    112.075 1951   63.221
## 1952         98.1 346.999      193.2        359.4    113.270 1952   63.639
# ggplot
g <- ggplot(longley, aes(x = Year, y = GNP)) +
  geom_line()+ # 折れ線
  geom_point()+ # マーカー
  labs(title = "1947~62年のアメリカのGNP")

# グラフ表示
g

# plotlyで表示
ggplotly(g)

9.3 箱ひげ図

箱ひげ図は五数要約の統計量（最小値、第一四分位数、中央値、第三四分位数、最大値）をプロットすることで分布の概形を表示する図である。

g <- ggplot(df, aes(x = class, y = japanese)) +
  geom_boxplot()+   # 箱ひげ図
  labs(title = "クラスごとの国語の点数の分布")


ggplotly(g)

9.4 バイオリンプロット

バイオリンプロットは、「カーネル密度推定」という手法で推定されたデータ分布の曲線を左右対称に描いたもので、異なるデータの分布を視覚的に比較できる。デフォルトでは最小値から最大値の範囲でプロットされ、geom_violin(trim = FALSE)にすると曲線の端はカットされない。

g <- ggplot(df, aes(x = class, y = japanese)) +
  geom_violin() + # バイオリンプロット
  labs(title = "クラスごとの国語の点数の分布")

ggplotly(g)

9.5 箱ひげ図とバイオリンプロットの併用

g <- ggplot(df, aes(x = class, y = japanese)) +
  geom_violin() +             # バイオリンプロット。これを先に描く。
  geom_boxplot(width = 0.3) + # 箱ひげ図を上書き。幅を0.3に指定
  labs(title = "クラスごとの国語の点数の分布")

g  # ggplotly(g)だと箱ひげ図が表示されない

10 層別化

ggplot()関数の中に入れるaes()関数のfillやcolorといった引数に質的変数を指定することで、色の塗り分けによる層別化ができる。

これにより３つ目の変数の情報を表現できる。

ggplot(df, aes(x = english, y = japanese, fill = class, color = class)) +
  geom_point() +
  labs(title = "英語と国語の得点の散布図をクラスごとに塗り分けたもの")

11 バブルチャート

散布図のサークルの大きさを、第三の変数に合わせて変更することもできる。

# ggplot
g <- ggplot(df, aes(x = english, y = japanese, size = math)) +
  geom_point(alpha = 0.5) + # alpha:バブルの透過度
  labs(title = "英語と国語の得点の散布図(バブルサイズは数学の点数を反映)") +
  scale_size(range = c(1,10)) # バブルサイズの範囲を設定

# 散布図表示
g

#plotly
ggplotly(g)

12 グラフの保存

図の保存にはggsaveを使用する

ggsave(plot = g, file = "bubble_chart.png")

（参考）事前に集計をしてからplotを描く場合

geom_bar()はstatという引数があり、デフォルトではgeom_bar(stat="count")と、観測値の数をカウントして表示するようになっている。そのためaes()にはxを指定するだけで自動的にy軸も描画された。

事前に集計を行って自分でx軸とy軸の変数を用意し、geom_bar(stat="identity")と指定することもできる。

# クラスごとに国語の平均点を集計
mean_ja <- df %>% group_by(class) %>% summarise(mean_ja=mean(japanese))
mean_ja
## # A tibble: 3 × 2
##   class mean_ja
##   <chr>   <dbl>
## 1 A        68.3
## 2 B        66.9
## 3 C        65.9
# プロット
mean_ja %>% ggplot(aes(x=`class`, y=mean_ja)) + geom_bar(stat="identity")

# pivot_longerで縦持ちデータに変え、科目ごとに平均得点を計算して描画する
g <- df %>%
  tidyr::pivot_longer(cols=c(japanese, english, math),
                      names_to="subject",
                      values_to="score") %>%
  dplyr::group_by(subject) %>%
  dplyr::summarise(mean_score=mean(score)) %>%
  ggplot(aes(x=subject, y=mean_score)) +
  geom_bar(stat="identity")

ggplotly(g)

13 参考文献

データ可視化については、様々な情報やテキストがある。また、グーグルの画像検索から自分が作りたいグラフ・イメージを探し、そのRコードを参考にするのも一つの方法である。