タグ：多重比較

2019年11月15日

【4-11】RでKruskal-Wallis（クラスカル・ウォリス）検定を行う方法

第4章は統計を扱います。

今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。

シロート統計学はEZRを使った統計分析をわかりやすく解説されています。

第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか？といった視点で進めていきます。

今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください！

今回はKruskal-Wallis（クラスカル・ウォリス）検定を紹介します

まずKruskal-Wallis（クラスカル・ウォリス）検定についてはハルさんのサイトをご参照ください。

また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。

1.準備

第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。

ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。

2.スクリプトファイルの作成

次にRのコードを書くためのスクリプトファイルを作ります。

3.データのダウンロード

今回もハルさんのサイトのデータを使わせていただきます。

デモデータ（Kraskal-Wallis検定）

この章ではRを使ってダウンロードしています。

download.file(url = “ファイルのURL”,
　　　　　　　 destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。

実際のコードは以下になります。
前回のコードのURL（" "の中）とdestfileのdata/以降を変更するだけでOKです。

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-kruskal-wallis-test.xlsx"
destfile = "data/demo-kruskal-wallis-test.xlsx"

download.file(url, destfile)

dataフォルダにダウンロードできたことを確認します。

4.データの読み込み

データを読み込みます。
今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。

View Fileでデータを確認します。

データが入っているセルを確認します。
B2からC92までデータが入っています（B2:C92と表記）
スクリーンショット 2019-11-12 7.51.15

次にImport Datasetでデータを取り込みます。

Import画面ではName, Sheet，Rangeを指定します。

Name：ハルさんのサイトではBWでしたがbw_kwとします（大文字・小文字は別物とされます）
Sheet：このExcelは1つしかデータがないのでDefaultのままでOK
Range：先ほど確認したB2:C92

Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。

#データの読み込み
library(readxl)
bw_kw <- read_excel("data/demo-kruskal-wallis-test.xlsx", 
                    range = "B2:C92")
View(bw_kw)

データが正しく入っていることを確認します。

これでデータの取り込みは完成です。

5.正規分布の確認

ハルさんのサイトに沿って正規分布の確認を行います。

正規分布の確認は【4-1】Rでt検定を行う方法で紹介しています。
以下のコードがイメージできない場合はまずこちらの記事をご参照ください。

まずヒストグラムで分布を確認します。
1つずつグラフを作るのは手間がかかるので、まとめてグラフを作ります。

ヒストグラムはgeom_histgram関数を使います。
グラフが重なっていると見にくいのでfacet.gridでカテゴリーごとにグラフを分けます。
ちなみに【4-10】Rで分散分析（一元配置分散分析）を行う方法で作ったヒストグラムのコピペでできてしまいます。変数名を変えるだけです。

library(tidyverse)
ggplot(data = bw_kw) +
  geom_histogram(aes(x = BW, fill = category), bins = 5) +
  facet_grid(category ~ .)

6.シャピロ・ウィルク検定
次にハルさんのサイトでは正規性の検定を行っています。
正規性の検討はshapriro.test関数を3つのカテゴリーそれぞれに行う必要があります。

今回はsplit関数とmap関数を使いまとめて行ってみます。
これも【4-10】Rで分散分析（一元配置分散分析）を行う方法で作ったコードと全く同じです。
変数名を変えるだけです。

コードを使うと他でも役に立つことがわかります。

bw_kw %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$BW))

まとめて結果が出ました。全てp < 0.05ですので正規分布であるという仮説が棄却されました。

7.グラフを作る
今回はノンパラメトリックの検定を行うので箱ひげ図を作成します。
箱ひげ図の作り方は【3-6】Rのggplot2で箱ひげ図を作るgeom_boxplot関数で紹介しています。

ggplot(data = bw_kw) +
  geom_boxplot(aes(x = category, y = BW))

もしEZRと同じようなグラフであればboxplot関数も使えるので併せて消化しします。

boxplot(数値 ~ グループ, data = ◯◯)
boxplot(数値 ~ グループ)

boxplot(BW ~ category, data = bw_kw)

8.Kruskal-Wallis（クラスカル・ウォリス）検定を行う

正規性が棄却されましたので今回はノンパラメトリックのKruskal-Wallis（クラスカル・ウォリス）検定を行います。

Kruskal-Wallis検定はEZRではkruskal.test関数を使います。

kruskal.test(目的変数 ~ factor(グループ), data = データ)

目的変数：BW
グループ：factor(category)
data = bw_kw

kruskal.test(BW ~ factor(category), data=bw_kw)

EZRと同じ関数を使っているので同じ結果になりました。

またMann-Whitney U 検定のときに紹介したcoinパッケージでも行うことができます。
もしcoinパッケージを使ったことがなければインストールします。

#coinパッケージのインストール
#1度でもしていればしなくてOK
install.packages("coin")

coinパッケージのkruskal_test関数はkruskal.test関数と同じ使い方です。
使う前にlibrary関数で呼び出します。
そしてグループにはfactor関数もしくはas.factor関数でfactor型に変更します。

kruskal_test(目的変数 ~ factor(グループ), data = データ)

library(coin)

kruskal_test(BW ~ factor(category), data=bw_kw)

Asymptotic（漸近的：正規近似）とあるので正確検定ではないようです。
そのためEZRと同じ結果になりました。
実は3以上のグループではなく2変数だと正確検定を行うことができます。

9.多重比較（Steel-Dwass法）を行う

EZRで使っていたSteel-Dwass法はEZR専用のパッケージなので、Rで行うにはNSM3パッケージのpSDCFlig関数を使う必要があります。

ちなみにEZRのSteel-Dwass法は正確検定ではなく正規近似を行います。
NSM3パッケージではデータ数が少なければ正確近似ができますが、データ数が大きくなると計算量が膨大になりすぎて実施できません。モンテカルロ法もしくは正規近似を行うことになります。今回の90例では正確近似はできませんでした。
正規近似と正確検定、Steel-Dwass法の使い方に関しては多重比較 Steel-Dwass 正規近似と正確検定： U 検定が基準に詳しく紹介されています。

#EZRと同じ正規近似で行う場合
pSDCFlig(目的変数, factor(グループ), method = "Asymptotic")

#正確検定を行いたい場合（データが多いと自動的にモンテカルロ法になる）
pSDCFlig(目的変数, factor(グループ))

目的変数：bw_kw$BW
グループ：factor(bw_kw$category)
＊グループはfactor型にする必要があるためfactor関数かas.factor関数を使います。
＊今回もdata = ◯◯ は使えないので $ を使います。

#Steel-Dwass法を行うためにNSM3パッケージをインストール
#1度でもインストールしてたらこの1行は必要なし
install.packages("NSM3")

#NSM3パッケージの呼び出し
library(NSM3)

#EZRと同じ正規近似で行う場合
pSDCFlig(bw_kw$BW, factor(bw_kw$category), method = "Asymptotic")

#正確検定を行いたい場合（データが多いと自動的にモンテカルロ法になる）
#結果に時間がかかるので注意！
pSDCFlig(bw_kw$BW, factor(bw_kw$category))

それぞれのp値は以下になりました。
大きな差はありません。
A-B：正規近似0.014、モンテカルロ0.0142
A-C：正規近似0.0255、モンテカルロ0.0249
B-C：正規近似0.9006、モンテカルロ0.902

10.ひとまずまとめ

今回はRでKruskal-Wallis（クラスカル・ウォリス）検定と多重比較であるSteel-Dwass法の紹介をしました。ハルさんのサイトではここまでですがBonferroniやHolm法でp値の調整をするにはどうしたら良いのか？やMann-Whitney U検定でさんざん問題になったタイと正確検定に関してはどうなったんだ？という方は続きがありますので興味がある方は続きをご覧ください。

次回は反復測定分散分析を紹介します。

11.他の多重比較の方法

一元配置分散分析でも紹介しましたが、他の多重比較にはBonferroniやHolm、Hochbergなどがあります。

一元配置分散分析の多重比較ではpairwise.t.test関数を使いましたが、Kruskal-Wallisではpairwise.wilcox.test関数を使います。簡単に言うと全ての組み合わせでMann-Whitney U検定を行いその結果を元にp値の調整を行います。
どの方法を使うかはp.adjust.method = で指定します。

pairwise.wilcox.test(目的変数, factor(グループ), p.adjust.method = "◯◯")

pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "bonferroni")
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "holm")
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "hochberg")

ただ、この方法だと警告が出ます。タイ（同順位）があるためです。
タイに関してはMann-Whitney U検定で紹介しました。

今回はずっと正規近似で行っているので正確近似で出すのであれば exact = FASEを追記します。
そうすればエラーは出ません。

pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "bonferroni", exact = FALSE)
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "holm", exact = FALSE)
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "hochberg", exact = FALSE)

11.どうしても正確検定を行いたい！！！という方は手計算が必要

多重比較は各組み合わせでMann-Whitney U検定を行いその結果を元にp値の調整を行います。
ただpairwise_wilcox_test関数が存在せず自分で手計算するしかありません。

そのためどうしても正確検定を行いたい場合はcoinパッケージのwilcox_test関数を使い、1つずつ組み合わせを行い、BonferroniやHolm、Hochbergそれぞれの方法で自分でp値を直接計算することになります。

この方法は統計ERさんの記事に詳しく紹介されています。
更にBonferroni、Holm、Hochbergって何をしてるの？という疑問も解決するかもしれません。
仕組みがわかればp値に2とか3とかかけるだけなので意外と単純です。

12.本当にまとめ
今回はRでKruskal-Wallis（クラスカル・ウォリス）検定と多重比較であるSteel-Dwass法の紹介をしました。ノンパラメトリック検定はタイと正確検定をするか正規近似を行うか？の問題がありますが、Mann-Whitney U検定のちゅんたろさんのツイートがここでも参考になると思います。

Chuntaro@Shuntarooo3
@MITTI12101 次に，タイがあってもある程度nがあれば正規近似でokです．nの数は文献により違いました．
『医学統計学ハンドブック』では2群の小さい方が8程度あれば実用上十分とあり，『医学への統計学』では2群あわせて30よ… https://t.co/8Av8FZk238
2019/11/01 10:40:35

そしてそもそもp値がより小さければより差が大きいわけではありませんので、正確検定か正規近似か？よりも研究デザインが適切かどうかの方が大切になるかもと言うのが個人的なイメージです。

次回は反復測定分散分析を紹介します！

13.今回使ったコード

今回使ったコードを紹介します。

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-kruskal-wallis-test.xlsx"
destfile = "data/demo-kruskal-wallis-test.xlsx"

download.file(url, destfile)

#データの読み込み
library(readxl)
bw_kw <- read_excel("data/demo-kruskal-wallis-test.xlsx", 
                    range = "B2:C92")
View(bw_kw) 

#正規性の確認
library(tidyverse)
ggplot(data = bw_kw) +
  geom_histogram(aes(x = BW, fill = category), bins = 5) +
  facet_grid(category ~ .)


bw_kw %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$BW))

#グラフの作成
ggplot(data = bw_kw) +
  geom_boxplot(aes(x = category, y = BW))

boxplot(BW ~ category, data = bw_kw)

#Kruskal-Wallis検定
kruskal.test(BW ~ factor(category), data = bw_kw)

#coinパッケージでKruskal-Wallis検定を行う場合
library(coin)
kruskal_test(BW ~ factor(category), data = bw_kw)


#Steel-Dwass法を行うためにNSM3パッケージをインストール
#1度でもインストールしてたらこの1行は必要なし
install.packages("NSM3")

#NSM3パッケージの呼び出し
library(NSM3)

#EZRと同じ正規近似で行う場合
pSDCFlig(bw_kw$BW, factor(bw_kw$category), method = "Asymptotic")

#正確検定を行いたい場合（データが多いと自動的にモンテカルロ法になる）
#結果に時間がかかるので注意！
pSDCFlig(bw_kw$BW, factor(bw_kw$category))


#他の多重比較

pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "bonferroni")
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "holm")
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "hochberg")

#正規近似を使えば警告は出ない。exact = FALSEを追加すれば良い
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "bonferroni", exact = FALSE)
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "holm", exact = FALSE)
pairwise.wilcox.test(bw_kw$BW, as.factor(bw_kw$category), p.adjust.method = "hochberg", exact = FALSE)

全記事の目次はこちらをご参照ください

2019年11月12日

【4-10】Rで分散分析（一元配置分散分析）を行う方法

第4章は統計を扱います。

今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。

シロート統計学はEZRを使った統計分析をわかりやすく解説されています。

第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか？といった視点で進めていきます。

今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください！

今回は一元配置分散分析を紹介します

まず一元配置分散分析（ANOVA）についてはハルさんのサイトをご参照ください。

また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。

1.準備

第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。

ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。

2.スクリプトファイルの作成

次にRのコードを書くためのスクリプトファイルを作ります。

3.データのダウンロード

今回もハルさんのサイトのデータを使わせていただきます。

デモデータ（ANOVA）

この章ではRを使ってダウンロードしています。

download.file(url = “ファイルのURL”,
　　　　　　　 destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。

実際のコードは以下になります。
前回のコードのURL（" "の中）とdestfileのdata/以降を変更するだけでOKです。

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-anova.xlsx"
destfile = "data/demo-anova.xlsx"

download.file(url, destfile)

dataフォルダにダウンロードできたことを確認します。

データが入っているセルを確認します。
B2からC22までデータが入っています（B2:C62と表記）
スクリーンショット 2019-11-08 14.26.26

次にImport Datasetでデータを取り込みます。

Import画面ではName, Sheet，Rangeを指定します。

Name：ハルさんのサイトではgripでしたがgrip_anovaとします（大文字・小文字は別物とされます）
Sheet：このExcelは1つしかデータがないのでDefaultのままでOK
Range：先ほど確認したB2:C62

Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。

library(readxl)
grip_anova <- read_excel("data/demo-anova.xlsx", 
                         range = "B2:C62")
View(grip_anova)

データが正しく入っていることを確認します。

これでデータの取り込みは完成です。

5.正規分布の確認

ハルさんのサイトに沿って正規分布の確認を行います。

正規分布の確認は【4-1】Rでt検定を行う方法で紹介しています。
以下のコードがイメージできない場合はまずこちらの記事をご参照ください。

まずヒストグラムで分布を確認します。
1つずつグラフを作るのは手間がかかるので、まとめてグラフを作ります。

ヒストグラムはgeom_histgram関数を使います。
グラフが重なっていると見にくいのでfacet.gridでカテゴリーごとにグラフを分けます。
facet.gridに関しては【3-4】Rのggplot2でヒストグラムを作るgeom_histogram関数で紹介しています。

library(tidyverse)
ggplot(data = grip_anova) +
  geom_histogram(aes(x = grip, fill = category), bins = 5) +
  facet_grid(category ~ .)

6.シャピロ・ウィルク検定
次にハルさんのサイトでは正規性の検定を行っています。
正規性の検討はshapriro.test関数を3つのカテゴリーそれぞれに行う必要があります。

今回はsplit関数とmap関数を使いまとめて行ってみます。
下の図は【4-1】Rでt検定を行う方法で使った方法ですが今回のコードと見比べてみてください。
コードがほとんど同じです。コードを使うと他でも役に立つことがわかります。

grip_anova %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$grip))

まとめて結果が出ました。すべて正規分布であることが棄却されませんでした。

7.等分散性の確認

次に等分散性の確認を行います。

3群以上で等分散性を検定するにはBartlett検定を使います。

bartlett.test(数値の列 ~ グループ)

data = を使えば$は必要ありません。どちらでも大丈夫です。

bartlett.test(grip_anova$grip ~ grip_anova$category)

bartlett.test(grip ~ category, data = grip_anova)

8.分散分析（ANOVA）を行う

それでは本題の一元配置分散分析と多重比較を行います。

実は一元配置分散分析と重回帰分析（t検定も）は同じ分析方法です（詳しく知りたい方は一般線形モデルで検索してみてください）。出てきた結果を分散分析の形式で表示するか重回帰分析の結果で表示するかの違いだけです。

そのため一元配置分散分析を行うには分散分析のaov関数か重回帰分析のlm関数で行います。
ただその後に行う多重比較ではlm関数が使えないものもあるので基本はaov関数で大丈夫です。

どちらの関数もまずモデルを作りサマリーを表示するというのが一般的な流れです。

全体的な流れは以下の通りです。

sumamry関数だけ結果の表示が変わるので注意が必要です。

aov関数の場合

モデル名：model_grip_anova（好きな名前でOK）

従属変数：grip

独立変数：category

データ名：grip_anova

model_grip_anova <- aov(grip ~ category, data = grip_anova)
summary(model_grip_anova)

lm関数の場合

モデル名：model_grip_lm（好きな名前でOK）

従属変数：grip

目的変数：category

データ名：grip_anova

model_grip_lm <- lm(grip ~ category, data = grip_anova)
summary.aov(model_grip_lm)

結果はどちらもEZRと同じになりました。

9.多重比較を行う

ハルさんのサイトでは多重比較でTukeyを紹介していますが他も紹介します。
Tukeyだけ関数が違います。

Tukeyの場合
Tukeyの場合は先程のaov関数の結果を使います。
lm関数の結果は使えませんので注意が必要です。

TukeyHSD(model_grip_anova)

Bonferroniの場合
pairwise.t.test関数はモデルの結果を使いません。

pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "bonferroni")

Holmの場合

pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "holm")

Hochbergの場合

pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "hochberg")

10.まとめ

今回は一元配置分散分析を紹介しました。
今回は割愛しましたが、一元配置分散分析を行うときはlongデータである必要があります。

Excelでデータ収集するときに注意するか、wideデータの場合はlongデータに直してから分析を行ってください。

longデータとwideデータについてはで【2-5】Rでデータを集計するのに便利なlongデータとgather関数で
紹介しています。

次回はKruskal-Wallis検定を行います。

（追記）グラフの作成
EZRでは検定と同時にグラフを作成しますが、Rでは直接作成します。
【4-1】t検定でグラフを使ったときのコードがそのまま使えます。

グラフを作るために必要な要約のコードもそのまま使います。

#データの要約
grip_anova_summary <- 
  grip_anova %>% 
  group_by(category) %>% 
  summarize(平均 = mean(grip),
              標準偏差 = sd(grip),
              '0%' = quantile(grip, 0),
              '25%' = quantile(grip, 0.25),
              '50%' = quantile(grip, 0.5),
              '75%' = quantile(grip, 0.75),
              '100%' = quantile(grip, 1),
              n = n())
grip_anova_summary

#グラフ化
ggplot(data = grip_anova_summary)+
  geom_bar(aes(x = category, y = 平均), color = "black", fill = "gray", stat = "identity") +
  geom_errorbar(aes(x = category, ymin = 平均 - 標準偏差, ymax = 平均 + 標準偏差), width = 0.1) + 
  labs(y = "grip") +
  theme_classic(base_size = 15) +
  scale_y_continuous(expand = c(0,0), limits = c(0,50))

コードの違いはわかったでしょうか？
t検定で使ったデータであるgripを今回のデータであるgrip_anovaに変更しただけです。
それ以外のコードの中身は全く一緒です！

11.今回のコード
今回使ったコードのまとめになります。

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-anova.xlsx"
destfile = "data/demo-anova.xlsx"

download.file(url, destfile)

library(readxl)
grip_anova <- read_excel("data/demo-anova.xlsx", 
                         range = "B2:C62")
View(grip_anova)


#正規分布の確認
library(tidyverse)
ggplot(data = grip_anova) +
  geom_histogram(aes(x = grip, fill = category), bins = 5) +
  facet_grid(category ~ .)


grip_anova %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$grip))

#Bartlett検定

bartlett.test(grip_anova$grip ~ grip_anova$category)

bartlett.test(grip ~ category, data = grip_anova)

#一元配置分散分析
#aov関数の場合
model_grip_anova <- aov(grip ~ category, data = grip_anova)
summary(model_grip_anova)

#lm関数の場合
model_grip_lm <- lm(grip ~ category, data = grip_anova)
summary.aov(model_grip_lm)

#多重比較
#Tukey
TukeyHSD(model_grip_lm)

#Bonferroni
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "bonferroni")

#Holm
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "holm")

#Hochberg
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "hochberg")

#データの要約
grip_anova_summary <- 
  grip_anova %>% 
  group_by(category) %>% 
  summarize(平均 = mean(grip),
              標準偏差 = sd(grip),
              '0%' = quantile(grip, 0),
              '25%' = quantile(grip, 0.25),
              '50%' = quantile(grip, 0.5),
              '75%' = quantile(grip, 0.75),
              '100%' = quantile(grip, 1),
              n = n())
grip_anova_summary

#グラフ化
ggplot(data = grip_anova_summary)+
  geom_bar(aes(x = category, y = 平均), color = "black", fill = "gray", stat = "identity") +
  geom_errorbar(aes(x = category, ymin = 平均 - 標準偏差, ymax = 平均 + 標準偏差), width = 0.1) + 
  labs(y = "grip") +
  theme_classic(base_size = 15) +
  scale_y_continuous(expand = c(0,0), limits = c(0,50))

全記事の目次はこちらをご参照ください

↑このページのトップヘ