多元素データの生成

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, 微量元素