初めに
近年の生物系の論文のグラフは生データを点で示し、検定結果とともに表示することが多くなっています。
このようなグラフはエクセルで作成するのは大変なのでRで作成したい、しかし記述量はできるだけ減らしたいというときは参考にしてください。
Rで最低限のグラフを出力し、パワーポイント等で整形する方法です。
ダイナマイト・プランジャー・プロット (平均値の棒グラフ＋エラーバー）は非推奨
コメントでご指摘いただきました。初版ではダイナマイト・プランジャー・プロットを使い記事を作成していましたが、このグラフは複数の測定値を集計して表示する際には適していない点があります。ダイナマイトプロットは以下のような図です。ダイナマイトプランジャーは、右の図のようなダイナマイトを発破する装置。
問題点としては、
・複数のデータを要約するグラフは、データが分布している領域を塗りつぶすべき。
・エラーバーは標準誤差や標準偏差、95％信頼区間などが使われるが、何が使われているのかわからず誤解を生む。正規分布をイメージしてしまうが実際わからない。
・平均値＋エラーバーでは情報を圧縮しすぎており、データの分布がわからない。
・元データを示しても平均値というただ一つの統計値が目立ちすぎ
参考
そこで今回はダイナマイトプロットの代替で書きます。
エクセルで読み込むファイルの作成
以下のような2列の表を作成します。
A1にはname, B1にはvalueと記述し、2行目以降のA列はサンプル名、B列は測定値を入れます。1列目をnameとvalueにしておけばこのページのコードをそのまま使えます。
今回はsample.xlsxというファイルのSheet1に入力しました。たくさん図を作りたい場合は、別シートで似たように作成してください。
サンプル名の並び順がそのままグラフの並び順になります。
| name | value |
|---|---|
| WT_Cont | 101.1 |
| WT_Cont | 80.7 |
| WT_Cont | 92.9 |
| WT_Cont | 106.8 |
| Mutant_Cont | 78.7 |
| Mutant_Cont | 73.7 |
| Mutant_Cont | 82.4 |
| Mutant_Cont | 76.1 |
| WT_Treat | 124.4 |
| WT_Treat | 129.5 |
| WT_Treat | 119.1 |
| WT_Treat | 125.7 |
| Mutant_Treat | 135.0 |
| Mutant_Treat | 121.6 |
| Mutant_Treat | 144.8 |
| Mutant_Treat | 120.0 |
作業フォルダーの変更。
読み込みたい エクセルファイルの存在する位置 に設定します。
Windowsの場合はRのウインドウの左上の ファイル ＞ ディレクトリの変更 です。Macの場合は command+d です。出てきたフォルダ選択の画面で作成したエクセルファイルの存在するフォルダーを選んでください。この場合は以下のsetwd()コマンドはいらないです。
R studioやJupyterLab、Linuxを使っている場合はエクセルファイルのあるフォルダーを右クリックして、パスのコピー、貼り付けて以下のかっこの中を書き換えます。ダブルクォーテーションを忘れずに。
Windowsの場合はバックスラッシュまたは円マーク \ を スラッシュ / に書き換えます。
setwd("C:/document/R")
ライブラリーのインストール
インストールされていないライブラリーのみインストールします。
以下の記事を参考にしました。
初回のみ以下の枠内を一気に貼り付けて実行します。
なんか色々出てくると思うので、終わるまで待ちます。結構長いです
">"が出てきたら終わりです
インストールに失敗したというエラーが出た場合は最新のRを入れなおしてください。
libs <- c("multcompView", "tidyverse", "ggpubr", "openxlsx", "svglite","export")
for(lib in libs) {
if(!require(lib, character.only = T)) install.packages(lib)
}
ライブラリーの読み込み
multcompViewはTukey検定で使い。tidyversはdataframeの操作に必要です。ggpubrはggplot2を使って便利に作図できます。openxlsxはエクセルファイルを開く際に使います。exportはパワーポイントを出力するのに使います。
Macではエラーになることがあります。OpenGLとrglに関連したエラーが出たら下のコマンドを試してみてください。
参考ページ:
library(multcompView)
library(tidyverse)
library(ggpubr)
library(openxlsx)
library(export)
options(rgl.useNULL = TRUE)
library(rgl)
library(export)
ファイルの読み込み
作成したエクセルを読み込みます。sample.xlsxというファイル名でSheet1というシート名から読み込む場合は以下のコードです。
適せん自分の用意したファイル名とシート名に変更してください。
data1 = read.xlsx("sample.xlsx", sheet = "Sheet1")
data1
Rは並び順を勝手に変えるので、固定します。
サンプル名を抽出し、並び順を指定するための配列を作成。
data1[ ,1]で表の1列目を指定しています。
order1 = unique(unlist(data1[ ,1]))
data1[ ,1] = factor(data1[ ,1], level = order1)
order1
グラフのy軸の範囲や検定結果を表示するためにデータの最大値を取り出します。
ymax = max(data1[ ,2])
グラフの作成
ggpubrパッケージでグラフを作ります。
参考ページ:
サンプルデータのような小規模なデータ（n < 10）では散布図と平均値がよいでしょう。
3行目はy軸の範囲を指定しています。０始まりにする理由がないときは３行目を消します。
1行目のx = "name", y = "value"は作成した表の列名に書き換えます。
fig = ggstripchart(data1, x = "name", y = "value", size = 3,
add = c("mean") , add.params = list(color = "red") )
fig = ggpar(fig, ylim = c(0, ymax * 1.1)) + scale_y_continuous(expand = c(0, 0, 0.1, 0))
fig
比較用として棒グラフ。できるだけ他の図を使います。add=c()の"mean"を"mean_sd"や"mean_se"にすると上で紹介したダイナマイトプロットになってしまいます。
fig = ggbarplot(data1, x = "name", y = "value", add = c("mean","jitter"))
fig = scale_y_continuous(expand = c(0, 0, 0.1, 0))
fig
箱ひげ図、サンプルサイズ n > 10 程度。あまり少ないデータだと四分位数の意味をなさないのと、一つのデータによってグラフの概形が大きく変わってしまいます。"jitter"の代わりに"dotplot"を使えば、dot plotという度数分布の一種を追加します。値が丸められる代わりに整列でき、最頻値がわかりやすいです。
fig = ggboxplot(data1, x = "name", y = "value", add = c("jitter"))
fig = ggpar(fig, ylim = c(0, ymax * 1.1)) + scale_y_continuous(expand = c(0, 0, 0.1, 0))
fig
さらにn数が多ければバイオリンプロット。add の中身の "jitter"を消せば、散布図は消えます。
fig = ggviolin(data1, x = "name", y = "value", add = c("jitter", "boxplot"))
fig = ggpar(fig, ylim = c(0, ymax * 1.1)) + scale_y_continuous(expand = c(0, 0, 0.1, 0))
fig
検定結果を追加
t検定か tukey検定の結果を追加できます。
1. T検定とともに出力
一行目は、検定したい組み合わせを書く c("A", "B")の形式で書きます。 増やす場合はカンマでつなぎ c("A", "B"), c("C", "B"), c("A", "D")のようにします。
5行目はt検定。 var.equal = TRUE でstudentのt-tesになり、var.equal = FALSE でWelchのt-testになります。
comp = list(
c("WT_Cont", "Mutant_Cont"),
c("WT_Treat", "Mutant_Treat") )
fig = fig + stat_compare_means(
comparisons = comp,
method = "t.test",
method.args = list(var.equal = TRUE) )
fig
2. Tukey test ともに出力
1行目の(value ~ name)は作成した表の列名に書き換えます。nameとvalueにしておけばそのままでよいです
anova = aov(value ~ name, data = data1)
tukey = TukeyHSD(anova)
tukey_res = multcompLetters4(anova, tukey, threshold = 0.05)[[1]]$Letters
tukey_res = data.frame(tukey_res)
tukey_res$name = factor(row.names(tukey_res),levels = order1)
tukey_res = arrange(tukey_res, name) #並び替え
tukey_res$group1 = tukey_res$group2 = 1
tukey_res
| tukey_res | name | group1 | group2 | |
|---|---|---|---|---|
| <chr> | <fct> | <dbl> | <dbl> | |
| WT_Cont | b | WT_Cont | 1 | 1 |
| Mutant_Cont | b | Mutant_Cont | 1 | 1 |
| WT_Treat | a | WT_Treat | 1 | 1 |
| Mutant_Treat | a | Mutant_Treat | 1 | 1 |
fig = fig + stat_pvalue_manual(
tukey_res, x = "name", y.position = ymax, label = "tukey_res",
position = position_dodge(0.8), size = 8
)
fig
テキスト部分の編集
パワポで編集もできますが、ここで編集する方法もあります。
Font sizeの変更
xy.textは軸の数値のフォントサイズ、xyは軸ラベルのフォントサイズ
fig = fig + font("xy.text", size = 20) + font("xy", size = 20)
fig
x軸ラベルを傾ける
fig = ggpar(fig, x.text.angle = 45)
fig
パワーポイントで保存
ファイル名は適当なものに変更してください。
別窓で出ているグラフのウインドウの大きさを変更すると保存されるファイルにも図のサイズが反映されます。保存したいサイズにマウスで合わせてもよいです。後にパワーポイント上で複数枚並べたい時などは、サイズを指定して揃えられます。
どちらか一方を実行してください
graph2ppt(fig, "sample.pptx") #表示されているサイズで保存
graph2ppt(fig, "sample.pptx", width = 8, height = 8) #縦横の長さを指定して保存
エクセルファイルの存在するフォルダに出力されるはずです。
SVGで保存
上のコマンドが動かない時や、パワーポイント以外（Inkscapeやイラストレーター等）で編集したいときはSVGで保存します。
ファイル名は自由に変更してください。(拡張子が.svgならSVG画像。.tiffならTIFF画像。.png にすればPNG画像として保存できます）最初に選択したフォルダー内にSVGファイルはできています。
ggsave("sample.svg", fig) #表示されているサイズで保存
ggsave("sample.svg", fig, width = 8, height = 8) #縦横の長さを指定して保存
パワーポイントで編集可能です。ドラックアンドドロップで図をパワーポイントに貼り付け、図を一度左クリックし、"Ctr+Shift+g"でグループ化の解除（windowsの場合）をして、ポップアップがでたら[はい]を選択します。Macの場合はCommand+option+Shift+gです。
その後2回ほどグループ化を解除すれば、ばらばらになるので、グラフの幅やフォントサイズ、色を変更可能です。
パワーポイントで加工
生データの点は枠線を太くすれば大きくなります。透明な四角形が重なってきているので、消すと、下のオブジェクトがクリックできるようになります。
散布図の平均値を横棒にパワポと手作業で置き換えて色付けするとこうなります。
