平均値の比較

Rとggplot2を使って平均値を比較してみましょう.

使うデータは前回作成した67元素のシミュレーションデータです.⇒ データ

 

d <- read.csv("67ele-sim.csv")

library(ggplot2)

ggplot(
  d, # 使用するデータフレームを指定
  aes(x=factor(area), y=Li) #指定したデータフレームからxとyを指定
  ) +
  geom_boxplot() # Boxplotを描画すると指定

# 背景や軸の色などを設定する
ggplot(d, aes(x=factor(area), y=Li)) + 
  geom_boxplot() +
  xlab("Area") + ylab("Concentration (μg/L)") + #x,y軸名を設定
  theme_bw() + # 基本的なテーマを白黒にする
  theme(
    panel.grid.minor = element_line(linetype = "blank"), #補助線を無し 
    axis.text = element_text(colour = "black"), # 軸のテキストを黒にする
    legend.position="bottom", # 凡例は下に表示
    text=element_text(size=14), # テキストのサイズを14に
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5) #x軸のテキストの位置を調整
  )

# 何のデータか分からないのでラベルを追加する
ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_boxplot(aes(fill=factor(area))) +
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]), #ラベルの位置を表示するためのデータフレームを作成
    aes(x=x, y=y, label=label), # 使用するデータを指定
    hjust=0, # ラベルの左端を水平の位置に合わせる
    vjust=1 # ラベルの上端を垂直位置に合わせる
    ) +
  xlab("Area") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )
# boxplotに平均値を追加する
ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_boxplot() +
  stat_summary( # 平均値を追加する
    fun.y="mean", # 算術平均を指定
    geom="point", # 点をプロット
    shape=21, # 縁のある丸
    size=2.5, 
    fill="white"
    ) + 
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  xlab("Area") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )
# バイオリンプロット
ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_violin(adjust = 1.2) +
  geom_boxplot(
    width=0.15, #ボックスの幅
    fill="Black", #ボックスの塗りつぶし
    outlier.alpha = 0 #外れ値の透明度を0にして表示しないように
  ) +
  stat_summary(fun.y="mean", geom="point", shape=21, size=2.5, fill="white") + # 平均値を追加する
     xlab("Area") + ylab("Concentration (μg/L)") +
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )
# 蜂群図
library(ggbeeswarm) # {ggplot2}のgeom_dotplot()でも描画できるが,{ggbeeswarm}を使った方がy値の違いを表現できる
## Warning: package 'ggbeeswarm' was built under R version 3.3.3
ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_violin(fill="grey90", adjust = 1.5) +
  geom_boxplot(width=0.15, fill="White", outlier.alpha = 0) +
  geom_beeswarm(
    colour="black",
    priority="density",
    cex=1.5,
    size=1
    ) +
    geom_boxplot(width=0.15, fill="White", outlier.alpha = 0, alpha=0.5) +
  stat_summary(fun.y="mean", geom="point", shape=21, size=3, fill="white") + 
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  xlab("Area") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )

Excelでこれらよりももっと簡単な棒グラフを作成して,x軸を地域間から採取月差を確認したいと思っても,12組のxとyの範囲を設定する必要が出てかなり面倒です.

Rを使って作図するようになると,x軸を他のパラメーターに変更するのはとても簡単に実行できます.

# xのデータをmonthに変更するだけで簡単に月別のデータを表示できる
ggplot(d, aes(x=factor(month), y=Li)) + theme_bw() +
  geom_boxplot() +
  stat_summary( # 平均値を追加する
    fun.y="mean", # 算術平均を指定
    geom="point", # 点をプロット
    shape=21, # 縁のある丸
    size=2.5, 
    fill="white"
    ) + 
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  xlab("Month") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )

さらに,地域×採取月のようなパターンの比較も簡単に作図できます.

# 地点ごと,採取月ごとの比較を行う

ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_boxplot() +
  stat_summary(fun.y="mean", geom="point", shape=21, size=2.5, fill="white") + 
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  facet_wrap(~month) +
  xlab("Area") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )

最後に,Boxplotで地域間の比較を67元素で一度に比較してアウトプットしてみます.

その際に,イオン価数,イオン半径,イオン化傾向など周期表の配置は元素間の関係を把握するために非常に重要です.そこで,周期表を模した配置してみます.これは{geofacet}というパッケージを使うと非常に簡単に作図することができます.

{geofacet}でしようする周期表上のグリッド( Grid_Period_Tbl_Ele (-ARI).csv )をまず作成・インポートします.


library(tidyr)
library(dplyr)
library(scales)

d.ele.melt <- tidyr::gather(d, "area", "month", 4:70) # Longフォーマットに変形
colnames(d.ele.melt) <- c("ID", "area", "month", "name", "Conc") # 列名を変更,元素は"name"にしておくこと
d.ele.melt$area <- as.factor(d.ele.melt$area)
  
library(geofacet)
PTEGrid <- read.csv("Grid_Period_Tbl_Ele (-ARI).csv") # カスタムグリッドは'code'と'name'にする必要がある

ggplot(d.ele.melt, aes(x=area, y=Conc)) + theme_bw() +
  annotation_logticks(sides = "l") + # l で左だけ,trblだと上下左右にtickが入る
  scale_y_log10( # y軸を対数軸に
    breaks = 10^(-7:7),
    labels = trans_format("log10", math_format(10^.x)) # y軸の表記を10^x
  ) +
  geom_boxplot(aes(fill=area)) +
  stat_summary(fun.y="mean", geom="point", shape=21, size=2.5, fill="white") +
  facet_geo(facets = ~ name, grid = PTEGrid, scales="free_y") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )
## You provided a user-specified grid. If this is a generally-useful
##   grid, please consider submitting it to become a part of the
##   geofacet package. You can do this easily by calling:
##   grid_submit(__grid_df_name__)

以下に示すように,各元素の地域別濃度を周期表上に表示することができます.

カーソルを移動させると,右側のウインドウで各元素のデータを確認することができます.

2017年10月30日