ブログ

サイト管理人のブログです。

ブログ一覧

ブログ始めました

テンプレートにblogが入っていたので、blogの項目を作りました。気が向いたら研究について何かしら書こうかと思います。

2017年08月29日ブログのカテゴリー:blog, その他

多元素データの生成

Rを使ってデータの可視化・解析を行うにあたって,自分たちが良く解析するデータを練習に使いたいと常々思っていました.まだ,論文にしていないものは公開するわけにもいかず… 必ずしも,解析したいようなデータがそろっていない…ということで,まずは実際の分析データに似たデータを生成するところから始めたいと思います.

ここでは,島根県内の水試料中に含まれる67元素の分析データの代表値から,多元素データを生成してみます.

一般的に,環境中の化学物質濃度は対数正規分布していることが多いので,幾何平均値と幾何標準偏差を計算したファイルを用意しました.さらに相関係数は,対数変換したのちにPearsonの相関係数を計算しました.

それぞれのファイルは以下の通りです (現在は非公開)

  • 幾何平均値:
  • 幾何標準偏差:
  • 相関係数行列:

正規分布と対数正規分布

library(mvtnorm)
# {MASS}の()では負の相関を扱えない

m = 67 # number of variable
g = 5
n = 60 # sample size
N = g *n

d.cor <- read.csv("cor.csv") 
m.cor <- as.matrix(d.cor[,-1])

GM <- read.csv("GM.csv")
GM <- GM[,-1]

GSD <- read.csv("GSD.csv")
GSD <- GSD[,-1]

#
df <- rmvnorm(n=N, mean=rep(0, m), sigma=m.cor)
# もう一度標準化しておく
df <- scale(
  df, 
  center = TRUE, # 平均値を0にする
  scale = TRUE # 標準偏差を1にする
  )

# for文を使ってセル[i, j]にlog(GSD[,j])を乗じて,log(GM[,j])を足して
# logから元に戻すために10^ の計算をする
for (i in 1:N)  {
  for (j in 1:m) {
    df[i, j] <- 10 ^ ( df[i, j]* log10(GSD[1, j]) + log10(GM[1, j])  )
  }
}

# 変数名を元素名に変更
colnames(df) <- colnames(GM)
d.ele <- as.data.frame(df)

# 後々比較を行うことを考えて,5つの地域と月の列を追加
area <- rep(1:5, 60)
month <- rep(1:12, 25)
df.ele <- cbind(area, month, d.ele)

summary(df.ele)
# データが確認出来たらcsvファイルとしてアウトプット
write.csv(df.ele, "67ele-sim.csv")
# Boxplotの作成
# ggplotでの作図はコードも長くなり,とっつきにくいですがより高品質なグラフを作ることができます
library(ggplot2)
library(reshape2)
library(scales)

d.ele.hist <- melt(d.ele) # 一つの変数として扱うために,データを列方向にまとめる
## No id variables; using all as measure variables
ggplot(d.ele.hist, aes(x=variable, y=value)) + 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(outlier.alpha = 0.6) + # boxplotを追加
  stat_summary(fun.y="mean", geom="point", shape=21, size=3, fill="white") + # 平均値を追加する
  xlab("") + 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=45, vjust=1, hjust=1)
  )
# 目的のグラフが作れたらファイルを保存する
# とくにWindows環境は透明度を上手く扱えないので,alphaの値を調製した場合はpngファイルで保存する
ggsave("boxplot.png", width=40, height=25, unit="cm", dpi=600)

# {ggplot2}の場合は,ggsave()で画像を保存できる.

環境分析を行った時にまず確認すべき事の一つは,測定値の妥当性です. 認証物質を測定して測定方法に問題が無いか確認しますが,正解が分からない環境試料の場合には,計算間違い一つで重大な影響を与える可能性もあり,より注意を払う必要があります.

  • 濃度レベルがこれまでの知見と大きな違いが無いか?
  • 外れ値はどうか?多すぎたり,偏っていたりしないか?
  • Oddo-Harkins則のように,(地球)化学的な知見から外れていないか?

以上の確認をして,場合によってはスペクトル干渉の有無を再確認し,他の分析条件を採用すべきか検討すべきでしょう.

適切なデータがあれば,以下のようなBoxplotが作成できます.

2017年09月29日ブログのカテゴリー:R, blog, ggplot, 微量元素

平均値の比較

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日

シンタックスハイライト

データサイエンス関係のブログを拝見していると,普段Rstudioで表示されるようにコードの視認性が高くなるようなスタイルが設定されています.HPBではRmarkdownで生成したhtmlをそのまま埋め込むことができないので少し放置していました.

このあたりの知識は全くなかったので,よく参考にしているHPのソースコードを見ると○○syntax highlightのような記述が・・・

このキーワードを検索してみたところ以下のサイトが比較的簡単に導入できそうな例を紹介しています.

JavaScriptでシンタックスハイライトする入力エリアを提供するライブラリを探した

このサイトの内容を参考にして,Rmarkdownで生成したhtmlファイルのソースコードからjavascriptとcssファイルを作成して,多元素データの生成の前半部分のコードに適用してみました.

 

library(mvtnorm)
# {MASS}の()では負の相関を扱えない

m = 67 # number of variable
g = 5
n = 60 # sample size
N = g *n

d.cor <- read.csv("cor.csv") 
m.cor <- as.matrix(d.cor[,-1])

GM <- read.csv("GM.csv")
GM <- GM[,-1]

GSD <- read.csv("GSD.csv")
GSD <- GSD[,-1]

#
df <- rmvnorm(n=N, mean=rep(0, m), sigma=m.cor)
# もう一度標準化しておく
df <- scale(
  df, 
  center = TRUE, # 平均値を0にする
  scale = TRUE # 標準偏差を0にする
  )

# for文を使ってセル[i, j]にlog(GSD[,j])を乗じて,log(GM[,j])を足して
# logから元に戻すために10^ の計算をする
for (i in 1:N)  {
  for (j in 1:m) {
    df[i, j] <- 10 ^ ( df[i, j]* log10(GSD[1, j]) + log10(GM[1, j])  )
  }
}

# 変数名を元素名に変更
colnames(df) <- colnames(GM)
d.ele <- as.data.frame(df)

# 後々比較を行うことを考えて,5つの地域と月の列を追加
area <- rep(1:5, 60)
month <- rep(1:12, 25)
df.ele <- cbind(area, month, d.ele)

HPBだけではjavascriptやcssファイルの場所を指定・保存できていないので,FFFTPで別途jsフォルダとcssフォルダを作成してその中に保存しました.

これはHPBのプレビューではコードにカラーが付いていることを確認できないので,何度か試行錯誤しましたが,なんとか上手くいったようです.

これまでの内容もすべて更新しました.モヤモヤしていたのがスッキリしました.これでもう少しブログの更新頻度を上げたいところです.

2017年11月09日

機械学習1 Random forest

我々も研究に機械学習を取り入れ始めています.機械学習(あるいは昨今のAIと思われているもの)というとDeep Learningのイメージが強いですね.Deep LearningはA Neural Network Playgroundというサイトで体験することが出来ます.このサイトではパラメーターを変化させてその挙動を見て学ぶことが出来ます.よく参考にさせていただいている有名なデータサイエンティストであるTJO氏のブログでも紹介されていますし,Deep Learning以外の機械学習の手法との比較をしています.

今回はその記事をなぞった内容で,上記NN Playgroundの二重らせんデータに対して何回かに分けて機械学習を試してみたいと思います.まずは二重らせんデータ生成し,用いました.本来はデータの一部のみを学習用データとするべきでしょうが,今回は全データを学習用データとして,TJO氏のブログにあるようにグリッド内の全点をテストデータとしています.また,私はggplot脳なので,{ggplot2}でアウトプットしていきたいと思います.

学習用データと試験データの作成

n <- 500
x <- seq(0.02, 4, by = 4/n)

noise1 <- rnorm(n, mean=0, sd=1)
noise2 <- rnorm(n, mean=0, sd=1)
noise3 <- rnorm(n, mean=0, sd=1)
noise4 <- rnorm(n, mean=0, sd=1)

L1 <- x + 0.05 * noise1 +0.02
L2 <- x + 0.05 * noise2 +0.02

r <- 0.05

X1 <- (1+r*noise3)*L1 * cos(10*L1/pi)
Y1 <- (1+r*noise3)*L1 * sin(10*L1/pi)

X2 <- -(1+r*noise4)*L2 * cos(10*L2/pi)
Y2 <- -(1+r*noise4)*L2 * sin(10*L2/pi)

df <- rbind(cbind(0, X1, Y1),cbind(1, X2, Y2)) 
df <- as.data.frame(df)
colnames(df) <- c("type", "x", "y")
df$type <- as.factor(df$type)

write.csv(df, "double spiral.csv")

px <- seq(-4.7, 4.7, 0.02)
py <- seq(-4, 4, 0.02)
grid <- expand.grid(px, py)
names(grid) <- c("x", "y")

NN PlaygroudでいうNoiseを加えています.実際はcsvファイルにアウトプットして毎度読み込んでます.{ggplot2}でNN Playgroud風にプロットするとこんな感じになります.

また,テストデータとしてgridという0.02刻みのデータを作成しました.

library(ggplot2)
theme_set(
  theme_bw() +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=12)
  )
)

ggplot(df, aes(x=x, y=y, colour=type)) + 
  geom_point(alpha=0.5) +
  scale_color_manual(values = c("darkorange", "darkblue")) 

Random forestを試してみる

我々は環境分析の研究室なので,環境汚染の度合いとその影響,発生源と環境汚染など目的変数と説明変数の関係性を評価したいというのが解析の目的になることが多いです.

そのため,変数の重要度を評価できるRandom forestという方法をよく使用しています.Random forestについての説明はWikipedia等見てください.

他にもその発展形のXgboostも変数重要度が計算できるので今後取り入れたいと思っているところです.

d <- read.csv("double spiral.csv")
d <- d[,2:4]
d$type <- as.factor(d$type)

px <- seq(-4.7, 4.7, 0.02)
py <- seq(-4, 4, 0.02)
grid <- expand.grid(px, py)
names(grid) <- c("x", "y")
library(randomForest)
## Warning: package 'randomForest' was built under R version 3.3.3
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
fit.rf <- randomForest(type~., d)
set.seed(8800)
tuneRF(d[-1], d$type, doBest=T)
## mtry = 1  OOB error = 0.7% 
## Searching left ...
## Searching right ...
## mtry = 2     OOB error = 1% 
## -0.4285714 0.05
## 
## Call:
##  randomForest(x = x, y = y, mtry = res[which.min(res[, 2]), 1]) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 0.5%
## Confusion matrix:
##     0   1 class.error
## 0 496   4       0.008
## 1   1 499       0.002
fit.rf <- randomForest(type~., d, mtry=1)

pred.rf <- predict(fit.rf, newdata=grid)
summary(pred.rf)
##     0     1 
## 94210 94661
d.pred.rf <- as.data.frame(cbind(grid, as.numeric(pred.rf)))
colnames(d.pred.rf) <- c("x", "y", "pred")
d.pred.rf$pred <- d.pred.rf$pred - 1 # 0,1にする

predict()関数でテストデータであるgridをd.rfというモデルで予測しています.予測した結果を対応する"x","y"に対して"pred"という変数で保存しています.

 

library(ggplot2)
theme_set(
  theme_bw() +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=12)
  )
)

ggplot(NULL) + # 異なるデータ数のデータフレームを重ねる時は初めにggplot(NULL) としておく
  geom_raster(aes(x=d.pred.rf$x, y=d.pred.rf$y, fill=d.pred.rf$pred), alpha =0.15) +
  geom_point(aes(x=d$x, y=d$y, colour=d$type), alpha = 0.5) +
  scale_fill_gradientn(colours = c("darkorange", "grey95", "darkblue")) + 
  scale_color_manual(values = c("darkorange", "darkblue")) +
  theme(legend.position = "none") +
  xlab("x") + ylab("y")

Random forestは決定木をベースにした解析法なので,TJO氏も指摘しているようにカクカクとした決定境界になっていますね.しかも,縦と横の線が組み合わさった境界になっています.また,(x,y)=(0,-0.8) では複雑なギザギザな境界が形成されて,yが-2.8付近では横に伸びた青の境界が形成されており,明らかに過学習していてます.

解析に使用する場合には注意したいところです.

2018年03月14日
» 続きを読む